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NOMENCLATURE 

ai,bi        Square  of  the  roots  of  the  characteristic  equation 
in  the  absorber  and  moderator,  respectively 

A^Bjj^Cjl  Constants  determined  from  the  boundary  conditions 

D.^  Parameters  for  empirical  fit  of  total  flux 

E  Energy,  or  Least  squares  error 

f(E>Qr)  Angular  flux 

fj/nbc)  Spherical  harmonics  component  of  angular  flux 

i,J,k,l,m,n  Summation  indices 

■fy /*(£=)       Associated  spherical  harmonics  (see  p.  224, 
reference  24) 

Q('r,-Q-)  Neutron  source 

Qj>„(;c)  Spherical  harmonics  component  of  neutron  source 

R^m  Elements  of  solution  vectors 

E  Field  position  vector 

dr.  Differential  volume  element 

r  Distance  perpendicular  to  cylinder  axis 

BS  Spherical  harmonics  component  of  scattering  oross 

section 

T  Matrix  containing  coefficients  of  the  A1 ,  Bif  and 

Cjl  as  given  by  the  boundary  conditions 

v  Neutron  speed 

X  Column  matrix  containing  the  Aj_,  Blt  and  C* 

x^  Elements  of  X 

z  Distance  along  cylinder  axis 

t-*i)?~Qx.  Radial  relaxation  constants  in  absorber  and  mod- 

erator, respectively 


iv 


&  Angle  between  jq.  and  the  z-  axis 

£^  Relaxation  constant  in  z  direction 

Jf0  Average  cosine  of  scattering  angle 

£  Total  macroscopic  cross  section,  Z^  +  £5 

I*  Macroscopic  absorption  cross  section 

Z$  Macroscopic  scattering  cross  section 

Ztr  Macroscopic  transport  cross  section,  Z.  -JfcZs 

P  Angle  between  r  and  the  projection  of  jql  onto  a 
plane  perpendicular  to  the  z-  axis 

$0:)  Total  neutron  flux 

&  Unit  vector  in  the  direction  of  neutron  motion 

d£-  Differential  solid  angle 


1.0  INTRODUCTION 

In  the  microscopic  theory  of  thermal  neutron  chain  reactors 
one  of  the  important  quantities  which  must  be  determined  is  the 
thermal  utilization  (8,24).  The  thermal  utilization  depends  on 
the  ratio  of  the  number  of  neutrons  absorbed  in  the  moderator  to 
the  number  absorbed  in  the  fuel.   In  a  homogeneous  reactor  this 
ratio  is  Independent  of  the  total  flux,  $ ,  but  In  heterogeneous 
reactors  it  depends  on  the  flux  distribution  in  both  the  fuel 
and  moderator. 

Several  methods  are  used  in  the  determination  of  the  absorp- 
tion ratio  and  the  suitability  of  a  given  method  depends,  among 
other  things,  on  the  optical  thickness  of  the  absorber.   For  a 
thin  absorber,  a  "first  flight"  calculation  can  be  used  (14). 
The  neutron  blackness,  which  is  defined  as  the  probability  that 
a  neutron  Incident  upon  a  body  will  be  absorbed  by  it  (19),  of 
an  absorber  depends  on  the  angular  distribution  of  the  entrant 
neutrons  and  this  can  be  determined  only  by  an  exact  solution  of 
the  transport  equation.   If  the  entrant  distribution  can  be 
approximated  as  that  described  by  diffusion  theory  (2),  Stuart's 
blackness  chart  (18)  can  be  easily  used  to  calculate  the  ratio 
of  absorption  in  the  moderator  to  absorption  in  the  fuel.  Numer- 
ical methods  can  be  used  ('6),  but  this  is  time  consuming  and 
each  problem  must  be  worked  individually. 

A  simple  method  for  calculating  the  absorption  ratio  is  to 
calculate  the  total  flux  distribution  In  the  lattice  from  the 
diffusion  theory  unit  cell  model.  As  Is  well  known  (1,24),  the 


ratio  determined  by  this  method  is  lower  than  the  experimental 
ratio.   The  experimental  flux  shows  a  greater  depression,  both 
in  the  fuel  and  in  the  immediately  surrounding  moderator,  than 
does  the  simple  calculation.   The  use  of  the  P,  or  higher  approx- 
imations  to  the  one  speed  transport  equation  reduces  this  discrep- 
ancy, but  does  not  account  for  temperature  effects  (24).  Never- 
theless, the  simple  methods,  just  because  they  are  simple,  are 
widely  used.   Fictional  cross  sections  are  sometimes  used  to  give 
better  results  between  theory  and  experiment  (26),  but  it  is  not 
clear  that  cross  sections  adjusted  to  give  agreement  in  one  lat- 
tice will  yield  agreement  in  another  lattice. 

Another  "parameter"  which  could  be  adjusted  to  correct  for 
the  faults  of  the  simple  theory  is  the  dimension  of  the  absorbing 
medium.  By  experimentally  determining  an  equivalent  radius  for 
different  absorbers,  moderators,  and  physical  dimensions,  a  param- 
eter would  be  available  which  could  be  used  to  give  the  correct 
results  by  the  use  of  a  simple  calculation.   The  general  useful- 
ness of  such  a  parameter  depends  on  how  well  the  equivalent  radius 
could  be  predicted,  as  a  function  of  lattice  parameters,  from  the 
results  of  a  relatively  few  experiments. 

The  theory,  including  appropriate  boundary  conditions  and 
computer  programs,  and  the  experimental  feasibility  of  obtaining 
an  equivalent  radius  for  the  P-^  approximation  by  making  flux 
measurements  in  assemblies  having  exponential  z-dependence  is  the 
subject  matter  of  this  work.  A  study  of  this  problem  for  the  P, 
approximation  is  being  made  by  Porath  (15). 


2.0  THEORY 
2.1  The  Spherical  Harmonics  Approximations 

Extensive  treatment  of  the  Boltzmann  equation  can  be  found 
in  many  references  (4,14,24).   In  this  work  only  a  brief  discus- 
sion of  the  general  equation  is  presented;  however  particular 
attention  is  placed  on  the  P,  approximation  to  the  spherical  har- 
monics component  form  of  the  Boltzmann  equation  in  cylindrical 
geometry.   It  is  assumed  that  the  diffusing  medium  is  homogeneous 
and  isotropic.  Only  the  monoenergetic,  time-independent  model 
for  neutron  transport  is  studied. 

The  general  Boltzmann  integro-differential  equation  is1 

v^t =  -•&■&*(£&,*>*)  -  Zb&,Bte)i&,£>£,*) 

(1 ) 

where  f(rrQ.,E,t)drdfl:dE  is  the  number  of  neutrons  in  the  volume 
element  dr.  having  energies  between  E  and  E  +  dE  whose  directions 
of  motion  lie  in  the  solid  angle  d£  about  £-,  multiplied  by  the 
neutron  speed  v;  Is(fc,t,E*E',£^a]!£u)/  is  the  cross  section  for 
changing  the  neutron  energy  and  direction  E&  into  the  range  dE7, 
dg£  at  E  and  •£■' due  to  collisions;  and 

%&■£**,*)   =  S*&j£Zskt,&*>-Qil-*±)  (2) 

is  the  cross  section  for  any  type  of  scattering. 

The  first  term  on  the  right  side  of  Eq.(l)  represents  the 


1  The  coordinate  system,  spherical  harmonics,  and  much  of  the 
nomenclature  used  in  this  work  nre  the  same  as  those  of  Weinberc 
and  Wigner  (24),  and  Kofink  (12).  6 


losses  due  to  the  straight  ahead  motion  of  the  neutrons.   The 
second  term  on  the  right  side  accounts  for  losses  due  to  absorp- 
tion and  scattering  out  of  the  phase  space.  The  gain  of  neutrons 
is  represented  by  the  third  and  fourth  terms  on  the  right  side  of 
the  equation.  The  third  term  accounts  for  all  types  of  sources, 
and  the  fourth  term  accounts  for  the  gain  resulting  from  the 
scattering  of  neutrons  into  the  phase  space  from  any  other  region. 

Davison  (4),  and  others  (8,14,24),  point  out  that  the  assump- 
tion of  a  monoenergetic  thermal-neutron  group  can  be  Justified 
only  for  slightly  absorbing  media  in  regions  away  from  sources 
and  boundaries.   This  assumption  has  yielded  results  which  agree 
reasonably  well  with  experiment,  however,  even  when  the  above 
restrictions  do  not  apply  (24).   The  energy- independent  thermal 
neutron  group  is  used  throughout  the  rest  of  this  work. 

Using  the  restrictions  that  the  medium  is  homogeneous  and 
isotropic  and  that  the  system  is  monoenergetic  and  time  independ- 
ent, the  Boltzmann  equation  reduces  to 

-£■%?&&)  ~   iPfoj)  +JWfli',J8./.^te/-*2)  +<?&£)  -  o  (3) 

Only  cylindrical  geometry  is  considered,  and  the  coordinate 
system  is  the  same  as  that  of  Weinberg  and  Wlgner  (24),  and 
Kofink  (12).  These  coordinates  are:  the  z-  axis  of  the  cylinder; 
the  distance  r-  from  the  axis  of  the  cylinder;  the  angle  e  between 
■Or   and  the  z-  axis;  and  the  angle  p   between  r-  and  the  projection 
of  a  In  a  plane  perpendicular  to  the  z-  axis. 

In  this  system  of  coordinates  the  Boltzmann  equation  for 
neutron  transport  is  given  by  Eq.(4). 


-g/j»A r./»c a  ?fifc»>fttf>    .    sine  sin  ^  M,y>9t*)      e.bS€>Wt>x,&>p) 

r  dr  i        r    dp  7y 

*  ('4) 

-  tf(t>%,0,f)   +  Jd*'#r>s0^ter-£)  +  90:,*)  =  o 

Meghreblian  and  Holmes  (14)  show  that,  since  the  scattering 
medius  is,  by  assumption,  homogeneous  and  isotropic,  the  cross 
section  £s(-fik*y  can  be  a  function  only  of  the  angle  60  between 
xj,  and*.  This  angular  dependence  can  be  conveniently  represented 
as  a  series  of  spherical  harmonics  of  the  first  kind  (11). 

Z6(s£*-fp    =  gs,Pf(coi00)  (5) 

Using  the  orthogonality  relation 

Jda  6gfe)^)  =  *&fek*>k&>  *  13k  *"*-*'  (6) 

the  scattering  cross  section  is 

2*  = Jte  Z/-3-'+4r)   -    4ffS0  (Y) 

The  average  cosine  of  the  scattering  angle  is  given  by 

XT  =    corgi    z,  J*&'«*4  z,(4+±l   _  Jx        Z'£tr 

^°  SWXtC&4>  ^   "  *s         .        ('8) 

To  realize  the  full  advantage  of  decomposing  the  scattering 

cross  section  into  spherical  harmonics,  the  same  is  done  for  the 

angular  flux. 

Ms*)  ~£*ster>*>k>M  (9) 

The  associated  spherical  harmonics,  P^Cfl-).  are  the  same  as 
those  of  Weinberg  and  Wigner  (24).  The  moments  are  given  by 

GMn>)  =  7F/^to6^)         do) 

The  £   =  0  and  7-1  terms  of  Eq. (10)' have  a  simple  physical 
interpretation,  but  the  higher  order  terms  have  no  simple  phys- 
ical meaning.  The  total  neutron  flux  at  r  is  obtained  by  inte- 
grating Eq.(9)  over  all  directions*.  Using  the  orthogonality 
relation,  the  total  flux  is  given  by  Eq.(ll). 


$(r,j)   =    [faffati)    =  Wfar,})  (11) 

The  z-  component  of  the  total  neutron  current  is 

/Jr>&  ~  fazwefciS)   *  f^^Aj)  (12) 

and  the  r-  component  of  the  total  neutron  current  is 

fr  t%$  s  jd-S-w  anp  %,-&)  =  -*g  K,!^)  -?itJnfi  ■    (13) 

To  obtain  the  spherical  harmonics  component  form  of  the 
Boltzmann  equation,  the  source  term  is  expanded  in  terms  of  the 
associated  spherical  harmonics  and  the  addition  formula 

is  used  with  Eq.(5)  and  then  substituted  into  Eq.(4).  Eq.(9)  is 
substituted  into  Eq.(4)  and  tho  resulting  equation  is  multiplied 
by  Py„(4-)dA  and  integrated  over  all  st.     Using  the  recursion 
relations 

si*'  & M»  «  Ll(zA$^0'&  -  (U-O^  ftletyiiJho  (15) 

the  spherical  harmonics  component  form  of  the  Boltzmann  equation 
is  obtained. 

4*r**«W£w  -U*-WCftp  -4,E»^^ 

*  ^  [£  -  ^J £V^  -  E,m  ^  £,r,p   .  £ £  £<< p  ( 16) 


where 


In  the  PL  approximation  to  the  spherical  harmonics  component 


form  of  the  Boltzmann  equation,  It  is  assumed  thai  ¥#   =  0  for 
/  ?   L.  For  L  equal  to  1  this  assumption  leads  to  the  same  equa- 
tions as  diffusion  theory,  hut  with  a  different  diffusion 
coefficient. 

In  the  P,  approximation  the  total  neutron  flux  Is  described 
by  a  second  order  differential  equation,  whereas  in  the  P*  approx- 
imation the  total  flux  must  satisfy  a  fourth  order  differential 
equation.  The  P,  approximation  should  give  a  more  accurate 
description  of  the  total  flux  in  regions  rear  boundaries,  sources, 
and  strong  absorbers,  where  the  flux  Is  a  rapidly  varying  func- 
tion. 
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2.2  Solutions  for  Particular  Conditions 

The  spherical  harmonics  component  form  of  the  Boltzmann 
equation  in  cylindrical  geometry  was  solved  for  two  different 
cases.   Both  cases  consist  of  a  two  region  medium  which  is  rad- 
ially symmetric.   In  this  work  the  central  region  is  referred  to 
as  the  absorber,  and  the  surrounding  region  as  the  moderator.  A' 
condition  of  symmetry  which  must  be  satisfied  for  both  cases  is 
that  there  be  no  current  around  the  z-  axis. 

The  first  case,  which  in  this  work  is  called  the  z-depen- 
dent  case,  is  a  two  region  medium  having  sources  at  the  z  =  0 
plane  and  having  one  surface  of  the  moderator  exposed  to  a  vacuum. 
The  second  case  is  the  unit  cell  problem  with  an  isotropic  source 
in  the  moderator. 

In  a  non-multiplying  medium  with  sources  at  z  ■  0,  GFlasstone 
and  Edlund  (9)  show,  that  in  regions  away  from  sources  and  the 
finite  end  of  the  cylinder,  the  z-dependence  of  the  total  neutron 
flux  can  be  described  by  a  decaying  exponential.   For  this  reason, 
and  because  of  the  simplicity  of  the  boundary  conditions,  it  was 
assumed  that  the  medium  extends  to  infinity  in  the  positive  z- 
dlrection. 

Letting  yf   -    1  -  4VSf/Laf+£l*.l  (jp) 

and  using  the  assumption 

+  fts  W«>  Av^r)  +  C4  UfiJo,)  Kjar)]  tUy  ( 19  ) 

in  Eq.(l6),  the  solution  for  the  z  dependent  case  with  linear 
anisotropic  scattering  is  derived  in  Appendix  A  and  is  given 


here  for  the  T-z   approximation. 
In  the  absorber  (central  region): 

bir>fi  '- Z/M^l^i"*  (20) 

In  the  moderator  (surrounding  medium): 

t  lu&DLtM&rt  *  (-ifCi^Uu^le^  (21) 

Letting  j=  !+£(£  + gO  (22) 

the  roots  of  the  characteristic  equation  in  the  absorber  are 

A2-  +o£  *  4fjD  +  a-i*8%V^y*;*J  =  a* 

**.*.7  (23) 

For  i  /  1  the  roots  of  the  characteristic  equation  in  the 
moderator  (Vfc2)  are  obtained  from  the  above  relations  by  replacing 
absorber  values  of  yo   and  Vj,  with  moderator  values.  For  i  -  1 

?-  gf  -  /fjli-^l-  i**fc  H,/3S- ^J    =■  bx  ( 24 ) 

Although  there  are  8  roots  in  each  region,  only  the  4  pos- 
itive values  are  retained  because  A  must  be  positive  in  order  to 
have  a  decaying  exponential  and  the  Bessel  functions  of  positive 
argument  are  not  independent  of  the  Bessel  functions  of  negative 
argument. 

The  RfK(<0  are  given  in  Table  1.  Except  for  1  =  1,  the 
R/fcCg.*')  are  calculated  from  the  same  relations  by  replacing 
absorber  values  of  /,  and  yt   with  moderator  values  and  by  replac- 
ing *<■   with  §^  and  a±   with  blt 

Using  the  relations  <r  =  -IX,  f   =  «,  l* /mfc'2S  +   l)/4Tr]^  z   r^ 
for  if  m   odd,  and  a^R2/  +  l)/4tj4  «  R^  for  ^«i  even,  where 
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Table  1 


Elements  of  Solution  Vectors 
for  the  z-Dependent  Case 


f 

m 

.1    =    1.2 

*W^y 

R^C-U) 

0 

0 

1 

0 

0 

1 

0 

3*oVaj 

0 

1 

1 

1 

3Wa  ,V2 

0 

-A/c(4V2 

2 

0 

(3Ai-aJ)NJ 

1 

5/iX/a4 

2 

1 

ocjXNjV6 

-2X/^Je> 

-5(2Ai-a4)M//6 

2 

2 

o»-NjV§/2 

(#+7)/«?/6 

-5*iV»4l^ 

3 

0 

3X(5^-3aj)Nj/5 

A 

yi(5^-a4)/a4 

3 

1 

3V3^(5Xi-aJ)N  /10 

-(3/f-7)/2«,V3 

-X(l5A2'-lla4)M/2/3 

3 

2 

9A*Jn  /V30 

A(3A*--7)/«jy30 

-5>*i.(3X2-a4)/a4/T0 

3 

3 

3^JN./2V5 

(X^)/^*^ 

-V5  ^4X/2a4 

aJNJ  =  ('aj-3y»yA-7y.)/(a  -7)  M  =  y^a^ 

'nie   RA^(?i^  are  OD"tained  from  the  R/„(<*i)  "by  using  moderator 
values  of  V0    and  ylf  replacing  a-^  with  blf  replacing  <*x  by  £i. 
and  then  multiplying  by  (-l)m. 
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1  is  the  square  root  of  -1,  the  coefficients  in  Table  1  can  he 
comp?red  with  the  coefficients  in  the  work  performed  by  Tralli 
and  Agresta  (23).   In  particular;  the  R52(«i),  R2o(°0>  R^t0^), 
and  R-,-z(  d4)  should  be  compared. 

The  unit  cell  model  is  frequently  used  in  heterogeneous 
reactor  calculations.   Some  of  the  assumptions  which  are  called 
for  in  this  model  are;  the  Blowing  down  density  is  constant  in 
the  moderator  and  zero  in  the  fuel,  the  fuel  lattice  is  broken 
up  into  a  number  of  identical  unit  cells,  and  the  fuel  elements 
are  so  long  that  the  z-dependence  of  the  flux  can  be  neglected. 

Using  the  same  notation  as  was  used  in  the  z- dependent  case, 
and  assuming  that 

i/r)  r  ^  ^(til^Lr)  iCLSM(<*)Wiri  ('25) 

the  solution  for  the  unit  cell  problem  with  an  isotropic  source 
term  and  linear  anisotropic  scattering  is  derived  in  Appendix  A 
and  is  given  here  for  the  P,  approximation. 


In  the  absorber: 


In  the  moderator: 


ftUr)  -  ^/L^Oc^JnJkx^  (26) 


WO   =    $  &/&■>  R^S.'O  +  Mflfc  **fc«I*J  +  f9  Sh*mo  ('27) 

Letting  <p    1    +    fcf$+-g>i)  (22) 

the  roots  of  the  characteristic  equation  in  the  absorber  are 

<*l  s  ^sLL+d-m^/ssf^J  (28) 

<**  *  7 

The  ot±   are  calculated  using  absorber  values  of  Y,   and  ft 
and  the  q   are  oaloulated  from  the  same  equations  by  using 
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moderator  values  of  V0   and  yA.  As  with  the  z-dependent  case, 
only  the  positive  roots  need  be  retained. 

The  Rf^oLiY  are  given  in  Table  2.   The  R/M(^')'  are  calcu- 
lated from  the  same  relations  by  replacing  absorber  values  of 
Vrf  and  ^  with  moderator  values. 


Table  2 

Elements  of  Solution  Vectors 
for  the  Unit- Cell  Model 


*     m 

Ra(Kj) 

.1    =   1.2 

R/*(  '«* ) 

0     0 

1. 

0 

1     1 

3^/=tjV2 

0 

2     0 

5N/4 

1 

2      2 

-5NjV678 

l/YS 

3     1 

3*jNjV5/8 

«^i/2V3 

3     3 

-3^NjV5/8 

^3/2^" 

Nj  =  1  -  3^V<</ 
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2.3  Boundary  Conditions 

The  exact  boundary  condition  at  an  interface  between  media 

is  that  the  angular  flux  shall  be  continuous.  This  implies  that 

the  moments  shall  be  continuous.  The  exact  condition  at  a  free 

surface  is 

^Ofifl)  -  O  for  4Et  inward,  (29) 

r  on  the  free  surface. 

In  the  spherical  harmonics  component  form,  an  infinite  number 
of  terms  are  required  in  order  to  satisfy  either  of  these  con- 
ditions. 

Three  methods  (4)  are  generally  employed  in  approximating 
Eq.(29).  The  first  is  collocation  wherein  f(r,£)  is  made  equal 
to  zero  at  the  required  number  of  points  &*•  Next  is  Mark's 
boundary  condition  where  one  imagines  the  vacuum  to  be  a  medium 
with  zero  scattering  cross  section  and  requires  the  angular  flux 
to  be  continuous.  The  last  consists  in  employing  a  set,  Z«M(&-), 
which  is  orthogonal  to  the  angular  flux  in  the  region  for  Jh 
inward,  and  then  making  J: f (r,£-)Z^  (^r)da-  equal  to  zero  for  r  on 
the  free  surface  and  jQ*  inward.   The  set  of  Z^n.)   is  usually 
chosen  from  the  spherical  harmonics  of  odd  order  £t    since  this 
set  automatically  includes  the  physical  requirement  that  the  in- 
ward neutron  current  is  zero.  These  are  called  Marshak's  bound- 
ary conditions  and  can  be  expressed  as 

jffc&){JrJa)<lSL  =o  for  4t  inward,  r  on  the  free  ('30) 

surface,  all  m7  /  odd  and  £  L 

where  Lis  the  maximum  /  in  an  odd  order  approximation  to  the 

spherical  harmonics  component  form  of  the  Boltzmann  equation. 
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In  the  Pt  approximation,  for  L  less  than  about  5  or  7, 
Marshak's  boundary  conditions  lead  to  the  best  convergence  (4). 
These  pre  the  boundary  conditions  which  were  used  in  this  worlc. 

As  is  shown  in  Appendix  A  and  as  can  be  seen  from  the  solu- 
tions for  the  f^trjz),  for  radial  symmetry 

Wr«>> =  ^ii<«*>  (31) 

Using  this  stipulation,    the  angular  flux  for  the  z-dependent  case 
becomes 

+    flb(r,}){;(3C0s',0 -i)  -  £/fy)v£  cos <j> sine coso  +  ^(r^fLtoszfistrfe 
1      ?3o  (r^)  i(scois&  -  3  cos  e)    -  £,  ftp  ±Ji  cospsme  (sce^a  -  -0 
i      ^(y^z^cosz.^  yrfe cos e   -  ^/ipf  vGr  tossfis/ne 
For   large  R,    the  outer  boundary  of  the   system,    Eq.(30) 
becomes 

J  (f (*#>£>% J£&*ete4d>   -   O  for  all  m,   /odd  and  6.  L .  (33) 

Substituting  Eq.(32)    into  Eq.(33)Nand  performing  the  indi- 
cated integrations 


(32) 


for  P 


10 


for  P 
for  P 


for  P 


11 
30 
31 


KA*)    t£fc£/*»£    =o  (34-a) 

£/*.*)  t£iZhfc£  -  £&/*,>)  +£&&<*,))   =  0  (34-b) 

?4//,p+;£ft/*,p-£Vi&tep*£Htep  =  o  (34-d) 

?#£//,,)  t  J-^  £#,£>  --  o  (34-e) 

-#•«*)  +  i&ftp  +i*y**>  +  £^£/'.p  --o  (34-f) 

Davison  (4)    shows  that  the  number  of  conditions  which  need 
to  be  satisfied  at  a  free   surface  is  equal  to   the  number  of  even 
order  moments,   N(j?=even),    and  the  conditions  to  be   satisfied  at 
an  interface  between  media  is  2N(^»even) .      The  use  of  either 


for  P^2 
for  PT, 
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Marshak'B  boundary  conditions  or  the  requirement  that  the  angular 
flux  shall  be  continuous  introduces  N(/«odd)  -  N(/=even)N  excess 
conditions. 

Most  applications  of  the  spherical  harmonics  method  have 
been  to  cases  where  the  number  of  odd  order  moments  was  equal  to 
the  number  of  even  order  moments.  For  cases  in  which  the  number 
of  odd  order  moments  is  greater  than  the  number  of  even  order 
moments,  Davison  (4)  suggests,  by  means  of  an  intuitive  argument, 
that  it  should  be  better  to  make  the  moments  for  all  m  with  Jl  < 
L-l  continuous  first  and  then  to  make  the  "predominately  normal" 
moments  of  order  L"  continuous.  The  same  consideration  should 
apply  to  the  choice  of  the  P^(h-)  in  Eq.(33). 

The  above  problem  does  not  arise  in  the  unit  cell  problem. 
Since  the  z- dependence  is  neglected,  the  angular  flux  should  be 
symmetric  about  any  plane  perpendicular  to  the  z-  axis.   This 
means  that 

Ur}e,<p)    =  -P(r,7r-9t4,)  ('35) 

From  an  inspection  of  Eq.(32)  it  is  seen  that  this  condition  can 
be  satisfied  if  only  such  ffw(r)  are  retained  as  have  /  +  m   even. 

In  addition  to  the  requirement  that  the  angular  flux  shall 
be  continuous  at  the  interface,  the  usual  condition  used  in  the 
unit  cell  problem  is  that  the  angular  flux  shall  be  symmetric 
about  the  cell  boundary  FL.   That  is, 

$(t<.>o><t>)   =  ?(£c>&,7r-0>)  ('36 ) 

Since  only  the  t/m(r)   with  A  m    even  are  retained,  an  inspection 
of  Eq.(32)  shoves  that  this  condition  can  be  satisfied  if  it  is 
required  that  Eq.(37)  be  satisfied. 
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fi,60  *£/0  =  i&>  -  °  (37) 

It  is  noted  that  the  number  of  even  order  moments  is  equal 
to  the  number  of  odd  order  moments  so  that  the  number  of  bound- 
ary conditions  is  equal  to  the  number  of  unknowns. 
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3.0  DISCUSSION  AND  RESULTS 
3.1  Consideration  of  Boundary  Conditions 

In  order  to  apply  the  equations  -which  have  been  derived  for 
the  P^  approximation  to  the  z-dependent  case,  an  appropriate 
choice  of  boundary  conditions  must  be  made.  For  a  two  region 
problem  with  known  cross  sections  and  dimensions,  Eqn''s  (20),  (21), 
and  (23)  show  that  there  are  13  unknowns  which  must  be  determined; 
X   and  4  each  of  the  A^,    B^,    and  Cj.   It  was  shown  that  there  are 
a  total  of  17  available  boundary  conditions;  10  from  the  require- 
ment that  the  angular  flux  be  continuous  at  the  interface,  6  from 
Marshak's  boundary  conditions  at  the  free  surface,  and  1  from  the 
source  condition.   In  order  to  have  a  meaningful  solution,  one  of 
the  unknowns,  say  A-^,  must  be  determined  by  the  source  condition. 
This  leaves  12  unknowns  to  be  determined  from  16  equations.  Of 
these  16,  any  8  of  the  10  obtained  by  matching  moments-  at  the 
Interface  plus  any  4  of  the  6  obtained  from  Marshal's  conditions 
form  an  appropriate  set  of  12  equations  in  12  unknowns. 

Two  procedures  for  solving  the  set  of  16  equations  are  to 
select  the  "best"  12  equations  or  to  retain  all  16  equations  and 
minimize  the  error  in  some  manner.  Most  work  with  the  p,  approx- 
imation has  been  to  cases  where  the  number  of  available  equations 
is  equal  to  the  number  of  unknowns  (3,12),  and  as  far  as  is  known 
to  the  author,  no  extensive  effort  has  been  made  to  determine 
which  equations  should  be  used  when  there  are  excess  equations. 
As  mentioned  previously,  Davison  ('4)  gives  an  intuitive  argument 
for  the  appropriate  choice.   Trail i  and  Agresta  (23)'  treat  a 
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problem  similar  to  the  previously  described  unit  cell  model  except 
that  they  consider  a  finite  cylinder.  For  their  problem  there  are 
12  unknowns  and  14  equations.   They  retain  all  14  equations. 
Which  of  the  two  procedures  gives  the  better  approximation  to  the 
exact  solution  of  the  transport  equation  is  not  known,  however  it 
is  much  easier,  and  is  consistent  with  the  method  used  in  the  P, 
approximation,  to  choose  the  "best"  boundary  conditions. 

Any  appropriately  chosen  combination  of  the  12  homogeneous 
equations  can  be  solved  by  arranging  the  equations  so  that  the 
following  matrix  equation  applies. 

TX  »  0  C38) 

Where  X  is  a  column  matrix  containing  the  A^,  Bj  ,  and  Cif  and  T 
is  a  square  matrix  containing  the  coefficients  of  the  unknowns  as 
given  by  the  12  equations.   Iii  principal,  the  application  of 
Cramer's  Rule  and  the  reduction  of  the  determinant  yields  an 
explicit  equation  for  the  A  which  make  the  determinant  of  T  equal 
to  zero.   Because  of  the  oscillatory  nature  of  the  Jm('r)  and 
Y^r),  there  may  be  an  infinite  number  of  these  AK.  For  each 
value  of  \,    any  11  of  the  12  equations  can,  in  principal,  be  used 
to  determine  all  the  x^  in  terms  of  one  of  the  x^  say  A]_.  The 
complete  solution  is  then  the  sum  of  all  these  solutions,  and  for 
the  z-dependent  case  has  the  following  form: 
In  the  absorber: 

fi*ft  i>   =  ! !/*  W«f JJkftfxr)  £Kly  (39 ) 

In  the  moderator: 
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where  the  R^n,  <x  ,  and  q   are  determined  from  the  same  relations  as 
before,  and  the  a*,  as  well  as  the  a|,  bJ,  and  cj  in  terms  of  a£, 
are  determined  from  the  "best"  boundary  conditions.   The  only  un- 
knowns in  the  above  equations  are  the  A*   and  these  can  be  deter- 
mined from  the  source  condition,  or  equivalently  from  the  spec- 
ification of  the  flux  at  some  position. 

In  order  to  get  some  idea  as  to  which  boundary  conditions 
should  be  used,  the  problem  was,  at  first,  simplified  by  consid- 
ering a  one  region  medium  having  a  free  surface,  and  containing 

the  origin.   For  this  problem  the  complete  solution  is  Eq.(40)  if 

k  v 

the  C  are  set  equal  to  zero.  Knowing  Br  from  the  source  condi- 
tion, A n   and  the  remainder  of  the  By  are  then  determined  from 
Marshak's  boundary  conditions.   This  reduces  the  number  of  pos- 
sible combinations  of  boundary  conditions  to  15. 

Instead  of  considering  all  15  possibilities  however,  it  was 
seen,  by  reasoning  along  the  lines  suggested  by  Davison  (4),  that 
any  combination  which  includes  the  equation  resulting  from  the 
use  of  P^q  in  Marshak's  boundary  conditions,  Eq.(30),  should  give 
a  poor  approximation  to  the  requirement  that  the  Inward  angular 
flux  be  zero,  because  f   is  the  predominately  "tangential  moment" 
of  order  3.  This  supposition  was  checked  by  using  what  was  con- 
sidered to  be  one  of  the  better  combinations  using  the  P,Q.   This 
was  called  Case  1  and  used  the  equations  resulting  from  the  use 
of  p10»  pn»  p30»  and  p31  in  EQ»(33).   By  comparing  Fig.  1  with 
Fig.  8  it  was  seen  that  Case  1  gives  a  much  poorer  approximation 
to  the  desired  condition  than  does  the  P-^  approximation.  Thus, 
all  combinations  containing  Eq.(34-c)  were  tentatively  rejected. 
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Of  the  remaining  5  combinations,  one  does  not  contain  the  equa- 
tion resulting  from  the  use  of  P,,  .   This  combination  defeats  the 
purpose  of  Marshak's  boundary  conditions  because  it  no  longer 
contains  the  stipulation  that  the  inward  neutron  current  be  zero. 
Of  the  remaining  4  combinations,  the  results  obtained  using  the 
P10,  P  ,  P,.,,  and  P32  in  Eq.(33)  should  be  similar  to  the  re- 
sults obtained  by  the  use  of  the  P1Q,  P,.,  P,2»  and  P-^,  since 

the  f   and  f-,,  are  the  "most  normal  moments"  of  order  3. 
31      33 

In  view  of  the  above  considerations,  3  additional  combina- 
tions of  the  P/»mwere  used  in  Eq.('33V.   These  are;  Case  2,  the 
P10,  P13,  P51,  and  P^J  Case  3,  the  P10,  P11,  P32,  and  P-^;  Case 
4,  the  P-q,  P^lt  P,2,  and  P^.   In  all  4  cases  only  the  first 
harmonic  in  Eq.(4o)was  used,  and  for  every  calculation,  B^  was 
arbitrarily  set  equal  to  1000.0.   Two  different  media  were 
studied,  iron  which  has  a  reasonably  large  absorption  cross  sec- 
tion, and  graphite  which  has  a  small  absorption  cross  section. 
Radial  dimensions  of  25.0  and  35.81  mean  free  paths, ££,  were 
used.   The  cross  sections  used  are  listed  in  Table  3. 

Table  3 

Cross  Sections  Used  in 
Study  of  Boundary  Conditions 

Medium Zi(c*-Q Me**)  ItthaOl 

Iron       0.222        0.933        1.144 
Graphite    2.83xl0~4      0.404        O.38I 

Some  of  the  results  are  shown  in  Figures  1  through  8.   The 
angular  flux  at  the  boundary  is  plotted  as  a  function  of  <f>   for 
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FIG.   3     CASE  2    (IRON,  2R  =  35.81) 
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FIG.  4     CASE2     (GRAPHITE,  2R  *  25.0) 
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FIG.  5     CASE  3     (GRAPHITE,  SIR" 35.81) 
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various  values  of  0.   In  addition,  the  figures  contain  a  curve 
for  which  the  ©-dependence  of  the  angular  flux  has  been  removed 
by  integration  over  the  ©-dependence  of  sl.      This  curve  is 

0  ( 4l) 

_vs  ;J§/4p  toy   -YF^^fy)  coi  3jz> 

where  R  is  the  radius  of  the  outer  boundary  and  Eq.(32)  has  been 
substituted  into  the  integrand.   Since  the  integral  of  sin©  from 
0  to  7T  is  equal  to  2,  Eq.(4l)  is  also  equal  to  twice  the  angular 
flux  averaged  over- the  ©-dependence  of_a.  Also  given  in  the  fig- 
ures is  the  partial  inward  flux 

J    J  Afc.j ,  e>t)  suede  ty     «  zvib  (A,%)  +  £  if  ft  (fy) 

+  *?«/«,*)-**  fi/J8tf)       ^2) 
The  parametric  dependence  of  f(R,  z,e,^>)  with  &   for  s  ;jr/t  is 
shown  only  in  Fig.  3  and  7,  the  cases  where  iron  was  the  medium 
under  consideration,  because  the  angular  flux  is  almost  symmetric 
about  e  -  V/l   for  the  graphite  medium.   There  was  a  slight  de- 
crease in  the  graphite  angular  flux  for  e  717/t.    ,  but  it  was 
small.   The  B±   and  the  value  of  X  which  describe  these  curves 
are  given  in  Table  4. 

These  curves  clearly  eliminated  Case  1  from  further  consid- 
eration and  possibly  indicated  that  Case  3  should  be  rejected 
also.   Cases  2  and  4  appear  to  satisfy  the  condition  that  the 
angular  flux  should  be  zero  for  Jl  inward  about  equally  well. 

In  order  to  proceed  further  it  was  necessary  to  consider 
the  two  region  problem  and  thus  to  make  a  choice  of  which  8  mom- 
ents should  be  matched  at  the  interface.   Since  Davison "s  (4) 
suggestion  (Case  2)  gave  good  results  at  the  free  surface,  his 
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Table  4 


Constants  Used  in  Figures  1  Through  8 


B±   =  1000. 

0 

Fig. 

B2 

*5 

B4 

X 

1 

-1.405xl0"25 

-9. 9 8 3x10" 20 

9.462xl0"56 

0.074597 

2 

-1.733x1c--27 

-1.302xl0"21 

-4.601xl0~4° 

0.079432 

3 

-2.586xl0"28 

-1.806xl0~18 

1.208x10-56 

0.70001 

4 

-3.776xlO"18 

-2.373xl0""14 

-1.337xl0"27 

0.10347 

5 

-3.363xl0-27 

-3.74lxl0~21 

-1. 840x10-57 

0.079375 

6 

-1.745xl0-27 

1.651xDO"21 

-1.049xl0"39 

0.079431 

7 

-5.546xl0-28 

3.262x10-18 

-9. 174x10" 37 

0.69998 

8 

0.0 

0.0 

0.0 

0.079498 

suggestion  for  which  moments  to  match  at  the  interface  was  heav- 
ily relied  upon.   As  previously  mentioned,  he  suggests  that  all 
moments  having  /  <  L-l  plus  the  predominately  normal  moments  of 
order  L  be  matched  at  the  interface. 

The  results  of  Case  2  indicate  that  the  moments  of  order  3 
which  should  be  considered  are  the  f^  and  f-,,.   Case  3  Indicates 
that  the  f,2  and  f-^  should  be  considered,  while  Case  4  suggests 
"the  fjlf  f-^,  and  f-^.  Neither  Case  2  nor  Case  3  conflicts  with 
Davison's  suggestion  concerning  the  6  moments  having  /  <  2,  but 
Case  4  does.   Since  the  spherical  harmonics  used  in  Case  4  all 
contain  at  least  some  normal  component,  all  "normal  moments"  were 
matched  first  and  then  the  remaining  moments,  starting  with  /  =  0 
until  the  required  number  of  equations  were  obtained.   Thus,  the 
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f20  and  f-^Q  were  not  made  equal  at  the  interface  for  Case  4. 

Using  these  rules  for  Cases  2,  3,  and  4,  the  total  neutron 
flux  distribution,  $  ,  was  calculated  for  a  central  region  of  iron 
and  a  surrounding  region  of  graphite.  The  cross  sections  of 
Table  3  were  used. 

The  distribution  calculated  for  Case  3  was  completely  unreal- 
istic, peaking  very  sharply  right  at  the  interface,  and  Case  3 
was  eliminated.  The  distribution  for  Cases  2  and  4  are  shown  in 
Fig.  9. 

It  was  seen  that  there  is  nothing  in  Fig.  9  which  indicates 
a  preference  between  Cases  2  and  4.   To  obtain  a  better  compar- 
ison, the  angular  flux  distribution  at  the  Interface  was  calcu- 
lated with  the  computer  program  described  in  Appendix  C.  A  rod 
radius  of  2.54  cm,  moderator  radius  of  50.0  cm,  and  the  cross 
sections  of  Table  3  were  used.  Only  one  harmonic  was  used  and 
A±   was  arbitrarly  set  equal  to  1.0.   The  results  are  listed  in 
Tables  5  and  6.  As  is  shown  in  Appendix  A,  f(r,z,a,*>)  = 
f(r,  z,e,-0>),  therefore  only  values  for  p   between  0  and  ir  are 
listed. 

This  calculation  clearly  showed  that  Case  2  gives  a  better 
approximation  to  the  requirement  that  the  angular  flux  be  contin- 
uous at  the  interface.   Therefore,  Case  2,  which  constitutes  the 
equations  suggested  by  Davison  (4),  were  used  In  the  remainder 
of  this  work. 
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Table  5 


Angular  Flux  at  Interface  for  Case  2 
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e 

P 

f (moderator) 

f (absorber) 

difference 

degrees 

degrees 

(n/cm  -sec/ 

(n/cm2- sec/ 

(n/cm2- sec/ 

steradian) 

stersdian) 

stersdian) 

0 

all  values 

2.923 

2.905 

0.018 

30 

0 

1.692 

1.741 

-0.049 

18 

1.756 

1.795 

-0.039 

36 

1.947 

1.958 

-0.011 

54 

2.251 

2.228 

0.023 

72 

2.635 

2.584 

0.051 

90 

3-040 

2.978 

0.062 

108 

3.400 

3.349 

0.051 

126 

3.666 

3.643 

0.023 

144 

3.824 

3.835 

-0.011 

162 

3.897 

3.936 

-0.039 

180 

3.917 

3.966 

-0.049 

60 

0 

1.387 

1.492 

-0.105 

18 

1.432 

1.518 

-0.086 

36 

1.602 

1.640 

-0.038 

54 

1.960 

1.938 

0.022 

72 

2.507 

2.438 

0.069 

90 

3.134 

3.046 

0.088 

108 

3.658 

3.588 

0.070 

126 

3.935 

3.914 

0.021 

144 

3.960 

3.998 

-0.038 

162 

3.860 

3.949 

-0.086 

180 

3.802 

3.906 

-0.104 

90 

0 

1.606 

1.606 

0.0 

18 

1.587 

1.587 

0.0 

36 

1.604 

1.604 

0.0 

54 

1.808 

1.808 

0.0 

72 

2.273 

2.273 

0.0 

90 

2.893 

2.893 

0.0 

108 

3.435 

3.435 

0.0 

126 

3.692 

3.692 

0.0 

144 

3.640 

3.640 

0.0 

162 

3.450 

3.450 

0.0 

180 

3.352 

3.352 

0.0 

Table  5   (continued) 
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& 

P 

f  ( moderator)1 

f  ('absorber) 

difference 

degrees 

degrees 

(n/cm2-sec/ 

(n/cm2-sec/ 

(n/cm2-sec/ 

sterodlan) 

steradlan) 

steradian) 

120 

0 

1.478 

1.373 

0.105 

18 

1.465 

1.379 

0.086 

36 

1.480 

1.442 

0.038 

54 

1.632 

1.654 

-0.022 

72 

1.986 

2.055 

-0.069 

90 

2.480 

2.568 

-0.088 

108 

2.959 

3.028 

-0.069 

126 

3.271 

3.292 

-0.021 

144 

3.373 

3.335 

0.038 

162 

3.346 

3.260 

0.086 

180 

3.318 

3.214 

0.104 

150 

0 

1.431 

1.382 

0.049 

18 

1.457 

1.418 

0.039 

36 

1.540 

1.529 

0.011 

54 

1.700 

1.723 

-0.023 

72 

1.937 

1.988 

-0.051 

90 

2.228 

2.290 

-0.062 

108 

2.525 

2.576 

-0.051 

126 

2.777 

2.800 

-0.023 

144 

2.953 

2.942 

0.011 

162 

3.051 

3.012 

0.039 

180 

3.082 

3.032 

0.050 

180 

all  values 

2.172 

2.190 

-0.018 

Table  6 


Angular  Flux  at  Interface  for  Case  4 
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9 

<P 

f (moderator) 

f (absorber) 

difference 

degrees 

degrees 

(n/cm^-sec/ 

(n/cm  -sec/ 

(n/cm  -sec/ 

steradlan) 

steradlan) 

steradian) 

0 

all  values 

2.357 

4.129 

-1.772 

30 

0 

0.667 

1.716 

-1.049 

18 

0.767 

1.816 

-1.049 

36 

1.062 

2.110 

-1.048 

54 

1.527 

2.575 

-1.048 

72 

2.115 

3.164 

-1.049 

90 

2.757 

3.805 

-1.048 

108 

3.368 

4.417 

-1.049 

126 

3.880 

4.928 

-1.048 

144 

4.249 

5.297 

-1.048 

162 

4.465 

5.514 

-1.049 

180 

4.536 

5.585 

-1.049 

60 

0 

1.618 

1.335 

0.283 

18 

1.675 

1.392 

0.283 

36 

1.870 

1.586 

0.284 

54 

2.243 

1.959 

0.284 

72 

2.785 

2.502 

0.283 

.  90 

3.403 

3.120 

0.283 

108 

3.947 

3.663 

0.284 

126 

4.295 

4.011 

0.284 

144 

4.427 

4.144 

0.283 

162 

4.427 

4.144 

0.283 

180 

4.410 

4.127 

0.283 

90 

0 

2.646 

1.858 

0.788 

18 

2.627 

1.839 

0.788 

36 

2.623 

1.836 

0.787 

54 

2.743 

1.955 

0.788 

72 

3.036 

2.248 

0.788 

90 

3.425 

2.638 

0.787 

108 

3.742 

2.955 

0.787 

126 

3.845 

3.057 

0.788 

144 

3.729 

2.942 

0.787 

162 

3.536 

2.748 

0.788 

180 

3.444 

2.657 

0.787 

Table  6  (continued) 


Z>t 


degrees   degrees 


f( moderator)    f( absorber)    difference 


(n/cm^-sec/ 

steradian) 


(n/cm2-sec/ 
steradian) 


(n/cm  -sec/ 
steradian) 


120 


150 


180 


0 

1.412 

18 

1.438 

36 

1.548 

54 

1.806 

72 

2.233 

90 

2.760 

108 

3.252 

126 

3.588 

144 

3.733 

162 

3.753 

180 

3.744 

0 

0.291 

18 

0.369 

36 

0.600 

.  54 

0.977 

72 

1.472 

90 

2.031 

108 

2.583 

126 

3.060 

144 

3.415 

162 

3.629 

180 

3.700 

all  values 

1.743 

1.302 
1.328 
1.438 
1.696 
2.122 
2.649 
3.141 
3.477 
3.623 
3.642 
3.634 

1.211 
1.289 
1.520 
1.897 
2.392 
2.951 
3.503 
3.980 
4.335 
4.549 
4.620 

3.120 


0.110 
0.110 
0.110 
0.110 
0.111 
0.111 
0.111 
0.111 
0.110 
0.111 
0.110 

-0.920 
-0.920 
-0.920 
-0.920 
-0.920 
-0.920 
-0.920 
-0.920 
-0.920 
-0.920 
-0.920 

-1.377 
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3.2  Practical  Solutions 

The  most  difficult  practical  problem  encountered  in  the  use 
of  the  P-^  approximation  to  the  Boltzmann  equation  for  neutron 
transport  in  the  z-dependent  case  is  the  solution  of  the  matrix 
equation,  Eq.(38).   In  the  z-dependent  case  the  determinant  of  T 
must  be  set  equal  to  zero  and  the  A*  determined  before  the  x£  can 
be  determined  in  terms  of,  say,  x£.  This  problem  is  not  encoun- 
tered in  unit  cell  calculations,  because  A  is  zero  and  the  equa- 
tions resulting  from  the  application  of  the  boundary  conditions 
are  Inhomogeneous.  Even  in  the  one  region  z-dependent  problem, 
however,  T  is  a  4x4  matrix,  each  term  of  which  is  a  complicated 
function  of  A ,  and  to  expand  the  determinant  to  obtain  the  char- 
acteristic equation  for  A  does  not  appear  profitable.  A  trial 
and  error  procedure  using  hand  calculation  methods  is  also  out 
of  the  question.   In  the  two  region  problem  16  different  Bessel 
functions,  some  having  as  many  as  10  different  arguments,  the 
various  R^,  plus  the  value  of  the  12x12  determinant  must  be 
determined  for  each  trial. 

The  computer  programs  described  in  Appendices  C  and  D  use  a 
linear  interpolation  procedure  to  select  succeeding  values  of  A 
in  an  attempt  to  make  the  determinant  of  T  equal  to  zero.   The 
"round  off  error"  in  the  calculation  of  a  12x12  determinant  can 
be  significant,  and  it  was  found  that  the  matrix  had  to  be  ordered 
so  that  each  diagonal  element  was  at  least  of  the  same  order  of 
magnitude  as  the  largest  element  in  the  same  column.   This  was 
particularly  important  for  the  equations  at  the  outer  boundary. 
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Because  the  dimensions  of  the  moderator  were  large,  the  equations 
at  the  outer  boundary  contained  both  very  large  and  very  small 
numbers-  (the  Im  and  Km  of  large  argument). 

Considerable  round  off  was  encountered  regardless  of  how  the 
matrix  was  ordered,  and  a  special  subroutine  was  developed  which 
expanded  the  numbers  to  18  digits,  performed  the  determinant  or 
matrix  calculations,,  and  then  reduced  the  answers  to  8  digits  for 
further  use  in  the  program.  This  eliminated  the  large  errors 
which  had  occurred  with  the  8  digit  arithmetic. 

It  was  not  expected  that  the  value  of  the  determinant  could 
be  made  exactly  equal  to  zero,  however  the  magnitude  of  the  error 
was  quite  unexpected.   The  results  of  a  typical  calculation  are 
shown  in  Table  7,  where  the  "slope"  is  defined  as  the-  rate  of 
change  of  the  determinant  of  T  with  A. 


Determinant 
of  T 


Table  7 
The  Determination  of  Lambda 

Slope        AxlO2 


3.210xl08° 
-3.394xl080 
-6.753xl078 

1.8l7xl077 
-1.228xl074 

5.629xl075 


-6.604xl082 
-6.473xl082 
-6.647xl082 
-6.653xl082 


-1.791x10 


83 


9.0000000 
10.000000 
9.4861780 

9.4757459 
9.4760192 

9.4760191 


It  is  seen  from  the  table  that  the  best  X  obtainable  still 
leaves  the  magnitude  of  the  determinant  quite  large.   The  slope 
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is  very  large,  however,  and  the  determinant  changes  sign  when  the 
'best1  X  is  changed  by  the  smallest  available  increment.  For  all 
practical  purposes,  then,  this  'best"  A  is  probably  the  correct  A. 

If  the  determinant  of  T  were  exactly  equal  to  zero  it  would 
make  no  difference  which  11  of  the  12  equations  were  used  to 
determine  the  x*  in  terms  of  x,  .  Since  the  determinant  was  not 
zero,  x-^  was  arbitrarly  set  equal  to  1.0,  each  of  the  12  equa- 
tions were  deleted  in  turn,  and  the  x*  were  calculated  from  the 
remaining  equations. 

The  maximum  difference  in  any  one  of  the  x,  was  approx- 
imately 1  in  the  fourth  digit,  and  this  occurred  in  only  one  case 
and  for  only  one  x..  The  maximum  difference  in  the  calculated 
total  flux,  j$ ,  at  any  position  was  approximately  5  in  the  seventh 
digit.  This  error  is  probably  less  than  the  error  due  to  round 
off  in  the  rest  of  the  program. 

Since  the  values  of  the  x.  are,  for  the  purposes  of  this 
work  (the  accuracy  in  $  is  more  important  here),  independent  of 
which  11  equations  are  used  to  determine  them,  the  'best1  value 
of  A  was  taken  to  be  the  correct  value. 
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3*3  Application  of  Theory- 
It  is  well  known  (1,8,24)  that  the  P-*  approximation  to  the 
Boltzmann  equation  for  neutron  transport  does  not  predict  as 
large  a  total  flux  depression  in  an  absorber  as  does  the  exact 
solution  to  the  transport  equation.  One  way  of  increasing  the 
depression  in  the  PL  approximation  is  to  use  a  fictional  absorber 
dimension  which  Is  larger  than  the  actual  dimension  by  an  amount 
Just  large  enough  to  give  the  same  ratio  of  average  moderator 
total  flux  to  absorber  total  flux  as  is  given  by  the  exact  solu- 
tion. What  is  more  Important,  from  the  practical  point  of  view, 
is  to  find  an  equivalent  radius  such  that  the  total  flux  described 
by  the  PL  approximation  describes  the  actual  distribution  found 
in  a  real  medium,  irregardless  of  whether  or  not  the  actual  dis- 
tribution is  described  by  the  exact  solution  of  the  transport 
equation. 

The  method  proposed  in  this  work  for  finding  this  equivalent 
radius  is  to  find  the  radius  which  makes  the  theoretical  expres- 
sion fit  the  experimental  points  as  accurately  as  possible.  A 
computer  program  which  is  described  in  Appendix  D  was  developed 
for  this  purpose.  The  fit  is  made  by  a  trial  and  error  procedure 
which  strives  to  minimize  the  sum  of  the  weighted  squares  of  the 
residuals  between  the  theoretical  expression  and  the  experimental 
values.   The  weighting  factor  is  the  reciprocal  of  the  standard 
deviation  squared  (22). 

In  order  to  test  this  program,  experimental  points  in  and 
near  an  iron  rod  embedded  in  a  graphite  medium  were  obtained  by  a 
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procedure  and  with  the  appartus  described  in  Appendix  B.  These 
points  were  not  intended  to  give  an  accurate  description  of  the 
neutron  flux  and  several  correction  factors  were  neglected. 
Foulke  (7)  gives  a  description  of  these  factors,  and  only  a  sum- 
mary of  some  of  these  effects  is  given  here. 

A'  hardening  of  the  neutron  spectrum  occurs  in  the  iron  rod 
due  to  preferential  absorption  of  the  low  energy  neutrons  in  the 
outer  regions  of  the  rod. 

Unequal  activation  on  different  sides  of  the  foils  can  occur 
since  foils  are  seldom  thin  for  neutrons  or  for  the  electrons 
produced  by  the  induced  radioactivity.  The  difference  in  activ- 
ity on  two  sides  of  a  foil  depends  on  the  flux  gradient  near  the 
foil.  As  can  be  seen  from  Figures  1  through  8,  this  effect 
should  be  more  pronounced  in  and  near  the  iron  rod  if  the  foil 
is  placed  perpendicular  to  the  z-   axis. 

A  foil  depresses  the  flux  in  a  region  around  itself  and 
since  the  foils  in  this  experiment  were  part  of  a  continuous 
strip  of  gold,  the  activity  induoed  at  any  point  in  the  strip 
was  influenced  by  the  presence  of  the  rest  of  the  strip. 

Since  gold  has  a  large  resonance  at  about  5  ev,  the  average 
energy  of  the  flux  which  induoes  the  activity  is  slightly  differ- 
ent than  the  average  energy  of  the  neutrons  present  in  the  system. 

In  addition  to  the  errors  introduced  by  the  neglect  of  these 
corrections,  other  errors  were  introduced  from  an  unknown  souroe. 
A  hump  in  the  neutron  flux  appeared  approximately  two  inches  from 
the  edge  of  the  rod  in  every  activation  and  the  cause  of  this 
hump  is  unknown. 
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The  results  obtained  by  fitting  the  theoretical  expression 
to  data  taken  from  Table  B-l  is  shown  in  Figures  10  through  15. 
The  effect  of  the  errors  mentioned  above  on  these  results  is 
probably  much  less  than  the  effect  caused  by  the  use  of  an  iron 
absorption  cross  section  which  was  much  too  large  (see  Appendix 
B).  Figures  10,  12,  and  14  show  the  best  fit  obtained  by  using 

1,  2,  and  3  harmonics,  along  with  the  equivalent  radii  which 
resulted.  Fig.  11  illustrates  the  relative  effects  of  using  1, 

2,  and  3  harmonics  for  the  actual  rod  radius  of  0.5  in.   Higher 
harmonics  would  presumably  give  better  results,  but  three  har- 
monics should  be  sufficient  if  the  measurements  are  not  made  too 
close  to  the  source.   Fig.  15,  which  uses  only  the  data  out  to 
the  hump  mentioned  above,  illustrates  the  type  of  fit  which  might 
be  expected  from  a  good  analysis  of  an  equivalent  radius.  The 
fact  that  the  fit  is  so  good  in  this  case  appears  to  be  the 
chance  result  of  a  lot  of  compensating  errors. 

It  is  to  be  noted  that  the  A  of  Eq.(39),  which  describe 
the  curves  in  the  above  mentioned  figures,  were  determined  solely 
on  the  basis  of  how  well  the  resulting  expression  for  the  total 
flux  fit  the  experimental  points.   They  therefore  presumably  have 
no  relation  to  the  A?  which  would  be  determined  by  a  source  con- 
dition. 

The  cross  sections  and  moderator  radius  which  were  used  to 
obtain  these  fits  are  discussed  in  Appendix  B.  No  extensive 
study  of  the  change  in  equivalent  rod  radius  as  a  function  of 
cross  sections  was  made.  However,  the  graphite  absorption  cross 
section  was  increased  from  2.53x10"^  cm-1  to  8.6x10"^  cm-1  and 
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the  resulting  increase  in  the  1  harmonic  equivalent  rod  radius 
was-  only  0.4$.  A  rough  estimate  of  the  effect  of  changes  in  the 
iron  absorption  cross  section  can  be  obtained  by  the  use  of  the 
P,  approximation  and  Fig.  13. 

The  ratio  of  the  flux  at  the  edge  of  the  rod  to  the  flux  at 
the  center  is  1.40  for  the  P-^  approximation  shown  in  Fig.  13. 
This  is  the  same  ratio  as  predicted  by  the  P^  approximation  for 
the  unit  cell  model  using  the  same  Xa,   of  0.259  cm"1.   The  exper- 
imental flux  ratio  in  Fig.  13  is  about  1.27.   Other  factors 
remaining  constant,  the  absorption  cross  section  would  have  to  be 
reduced  to  0.17  cm-1  in  order  for  the  P,  unit  cell  model  to  give 
the  Bame  ratio.  Even  if  the  actual  iron  absorption  cross  section 
(0.229  cm""l)  had  been  corrected  to  an  effective  neutron  temper- 
ature (14)  and  used  in  these  calculations  it  is  doubtful  if  the 
1  harmonic  equivalent  radius  would  have  been  as  large  as  the 
actual  rod  radius.   It  is  quite  possible,  however,  that  the  3 
harmonic  radius  would  have  been  at  least  as  large  as  the  true  rod 
radius. 
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3*4  Conclusions 

When  the  P^  approximation  to  the  Boltzmann  equation  for 
neutron  transport  is  used  in  cylindrical  geometry  which  does  not 
have  symmetry  about  planes  perpendicular  to  the  z-  axis,  appli- 
cation of  the  usual  boundary  conditions  gives  rise  to  excess 
equations  (4).  Of  the  combinations  studied  here,  the  desired 
conditions  were  most  nearly  met  by  matching  all  moments  for  /  <  3 
and  the  normal  moments  of  order  3  (the  f^  and  ty^)   at  the  inter- 
face along  with  the  use  of  the  P1Q,  P1;L,  P51,  and  P33  in 
Marshals' s  boundary  conditions. 

The  angular  flux  obtained  by  matching  all  moments  having 
m  f   0  plus  the  f00  and  f1Q  at  the  interface,  along  with  the  use 
of  the  P1;L,  P51,  P-^2,  and  P-j^  in  Marshak's  boundary  conditions 
gave  a  good  approximation  at  the  outer  boundary,  but  was  found 
lacking  at  the  interface.  Even  so,  the  total  flux  obtained  by 
using  these  conditions  shows  a  greater  depression  in  the  absorber 
than  the  flux  obtained  by  using  the  preferred  boundary  conditions, 
and  may  more  nearly  approximate  the  total  flux  as  described  by 
the  exact  solution  of  the  Boltzmann  equation. 

The  equivalent  radii  obtained  by  fitting  the  theoretical 
curve  to  the  experimental  data  were  in  error  because  of  the  neg- 
lect of  correction  factors,  the  use  of  inaccurate  cross  sections, 
and  uncertainty  in  the  "cleanness"  of  the  counted  activity.   In 
spite  of  these  errors  the  results  indicate  that  with  a  good  exper- 
imental determination  of  the  flux  and  with  the  correct  cross  sec- 
tions, a  good  determination  of  the  equivalent  radius  is  possible. 
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3.5  Suggestions  for  Further  Study- 
Much  of  the  experimental  work  in  the  determination  of  an 
equivalent  rod  radius  is  still  to  be  done.   In  addition  to  cor- 
recting the  faults  of  the  experimental  method  used  in  this  work, 
a  study  of  the  effect  of  changes  in  both  scattering  and  absorp- 
tion cross  sections  on  the  equivalent  radius  should  be  made. 
This  analysis  could  be  made  rather  easily  with  slight  modifica- 
tions of  the  computer  program  described  in  Appendix  D,  and  it 
would  show  which  cross  sections  should  be  known  with  the  most 
accuracy.   It  appears,  for  example,  that  the  exact  value  of  the 
moderator  absorption  cross  section  is  unimportant  as  long  as  it 
is  small  compared  to  the  scattering  cross  section. 

Since  the  higher  harmonics  die  out  with  increasing  distance 
from  the  source,  there  should  be  a  region  where  the  total  flux 
is  described  almost  entirely  by  the  first  harmonic.  By  making 
parallel  measurements  along  the  axis  of  the  rod,  one  should  be 
able  to  determine  whether  or  not  the  contribution  of  the  higher 
harmonics  is  significant.  If  the  flux  can  be  adequately  described 
by  the  first  harmonic,  a  conriderable  amount  of  computer  time  can 
be  saved. 

If  the  total  flux  can  be  described  with  the  use  of  only  the 
first  harmonic,  there  is  a  possibility  that  the  z-dependent  model 
can  be  replaced  by  an  equivalent  unit  cell  model.   The  require- 
ment at  the  cell  boundary  is  that  the  fn,  f^,  and  f^   be  zero. 
An  estimate  of  an  appropriate  cell  boundary  could  be  obtained  by 
plotting  the  fn,  f31,  and  f^   from  the  z-dependent  model  and 
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determining  the  position  at  which  this  requirement  is  most  nearly- 
satisfied.  If  the  unit  cell  model  can  he  used,  the  determination 
of  the  equivalent  rod  radius  would  be  greatly  simplified. 

Although  the  boundary  conditions  used  in  this  work  give  a 
good  approximation  to  the  desired  conditions,  one  further  pos- 
sible combination  might  be  considered.   Case  4  was  treated  rather 
inconsistently  in  that  some  of  the  moments  which  were  used  at  the 
outer  boundary  were  neglected  at  the  interface.   Since  Case  4  did 
give  a  good  approximation  at  the  outer  boundary,  better  results 
might  be  obtained  by  neglecting  the  f10  and  f,Q  at  the  interface 
as  well  as  at  the  free  surface. 
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(A-l) 


APPENDIX  A 

Solution  of  the  P3  Approximation  to  the 
Boltzmann  Equation  for  a  Z-Dependent 
Problem  and  for  the  Unit  Cell  Problem 

The  spherical  harmonics  component  form  of  the  Boltzmann 
equation  for  monoenergetic  neutrons  at  steady  state  in  a  homo- 
geneous, isotropic  media,  for  cylindrical  geometry,  using  the 
?fm(&)   and  nomenclature  of  Kofink  (12)',  is: 

-(I-   ifnlUty)    +  Q*n(*p     =  ° 
where 

*      — a^AiJ —  ■*"   — 2717=1) liW' 

fm        Wfl *7=i 

For  the  z-dependent  case  it  was  assumed  that  Q^,*  =  0  and 

that  the  neutron  sources  are  in  the  z  =  0  plane.  Noting  the 
following  relations  for  the  modified  Bessel  functions  of  integer 
order  m  (25) . 

Hr+  »&ImJ<*r)   =  t£r-"Tllm.xM    =  «l>r) 
[£+  aPlKmtl(4i)   =  E/r"  TJCjf-r)  =  -«&£>■)        U"3) 
it  is  seen  that  a  solution  of  the  form 

may  satisfy  the  homogeneous  part  of  Eq.(A-l). 

Setting  Qim  equal  to  zero,  placing  Eq.(A-4)  into  Eq.(A-l), 
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(A-6) 


(A-7) 


defining 

if    -  i  -  <WVi)X  (A-5) 

and  equating  coefficients  of  like  terms  yields  the  following 
relations. 

Comparison  of  the  above  relations  shows-  that 

Therefore,  only  the  relations  for  the  R/m(«*y  are  studied. 

If  there  is  no  net  flow  of  neutrons  around  the  axis  of  the 
cylinder,  the  condition  of  symmetry  to  be  satisfied  is 

?(*?,*>*)   *    $(r,z,6r<p)  (A-8) 

This  means  that  only  such  combinations  of  f?w('r,z)p,w(ja)s  as  con- 
tain cos(m^)  need  be  retained.  For  any  m  /  0,  the  angular  flux 
contains  pairs  of  the  form 

By  requiring  t^mMm   (-l)mf//Hf  the  symmetry  condition  is  satisfied. 
This  also  requires  that 

h-"     ~     (~$mR*.m  (A-10) 

Using  the  above  restrictions  for  the  V^   approximation  to  the 
Boltzmann  equation  (R,m  =  0  for  /  7   3),  and  assuming  linear  aniso- 
tropic scattering  (Bj  m   0  for  J  7   2),  the  following  set  of  equa- 
tions were  obtained  from  the  relations  for  the  R/M(a)   in  Eq.(A-6). 
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Using  a  Crout-type  reduction,  the  above  matrix  was  arranged 
so  that  the  first  four  rows  were  the  same  aB  above,  the  first 
four  columns  of  rows  5  through  10  contained  zeros,  and  the 
remainder  of  the  matrix  was  as  shown  below. 


10  0     0 

0     <f  o     0 

0     0  A     0 

0     0  0     0 

0     0  0     0 

0     0  0   -*£ 


XE/hYZ  DVS/2A 

(i<*-A*)E/AV§'        -VgX(^  +  D)/2A 
AE  D-2Xi 

0  C 

B  -A(8o^+9Az+I>.35) 
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R, 


00 


CA-12) 
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where 

A-  =  az+  £-7  D  =  <*£-3^>W>£     E  =  3>i+7     (A-13) 

B  =  (8>i  +  7)  (**+>*)  -35  n  (1-14) 

C  «  9(1*+ X*)2  -(35+27ye  >i  +  28Yo)(<*2+  ^)  +  l05yoyi       (A-15) 
Applying  Cramer "s  Rule,  and  setting  the  determinant  of  the 
matrix  equal  to  zero,  the  characteristic  equation  is 

ABC:  =  0  (1-16) 

Designating  the  roots  of  this  equation  in  the  following  manner, 
Az  +  c£    =^3[l-(i -iosy0 Xi/zsfj'*]  -  9£ 

aS  4  =  f  jEi+a-mUiM^y  =  ai 


^i 


X  +  wf  =  7 


(1-17) 


A"  +  «<J  tr  3^/^+7)  =  a4 
where  5  =  i  +  Y0  ( ^  +  ^sVt)  (1-18) 

the  general  solution  for  the  moments  is 

The  R^m,  as  given  in  Table  a-1,  were  obtained  by  separately 
placing  the  roots  into  the  matrix  and  arbitrarily  setting  one  of 
the  Rjn  equal  to  1. 

For  a  medium  as  shown  in  Pig.  A-l,  the  cylinder  is  unbounded 

in  the  positive  z-  direction  and  the  boundary  conditions  are 

Llm  fj»w(r,z)  =0    (all  r)  (A-20) 

z-»+«>  ' 

f,„('0,z)  is  finite  (1-21) 

In  addition,  the  moderator  «±   will  be  imaginary.  Desig- 
nating the  roots  of  the  characteristic  equation  for  the  moderator 
by  bi  and  using  ^   instead  of  ctit    the  b^  ^  ,  and  R^C^V  for 
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Table  A-l 

Elements 
for  the 

of  Solution  Vectors 
Z-Dependent  Case* 

i 

m 

.1  =  1.2 

Ra,(«>) 

R^(^) 

0 

0 

1 

0 

0 

1 

0 

3/oX/a 

0 

1 

1 

1 

3y«otj/ajV2" 

0 

-A/*«V2 

2 

0 

(3*-^ 

1 

5>lVa4 

2 

1 

otjX^  V6 

-2Vo<s/6 

-5(2Al-a4)M/V6 

2 

2 

^N^/2 

(Ai+7)/«<11V6 

-5nVa4V5" 

3 

0 

3A(5A2-3aJ)NJ/5 

A 

^i.(5A2-a4)/a4 

3 

1 

3V3ctj(5A2-aj)Nj/lO 

-fVf-Dfrxjfi 

-A(15AMla4)M/2V3 

3 

2 

9AotJNj/V30 

X(3Az-7)/c<tV3b 

-5y1(3Ai-a4)/a4V30 

3 

3 

34-NJ/2V5 

('X*+ 7)7204^ 

-V5  >ioC,A/2a4 

Vj  =  (ar3V<'yi"7Y,')/(aj"7)  M  =  *i/a4«4 

#The  Rjjt'fo)   are  obtained  from  the  R^f'o^)'  by  using  moderator 
values  of  a;  and  V±,  replacing  a-  with  b,  ,  replacing  *,  by  £t 
and  then  multiplying  by  (-l)m. 


59 


Absorber 


Vacuum 


Moderator 


Vacuum 


Fig.  A-l  Geometry  for  the  Z- Dependent  Case 


i  ?   1  are  calculated  In  the  same  manner  as  the  corresponding 
relations  in  the  absorber.  For  1=1  the  RfJUfc}   are  given  in 
Table  A-l  and  gL  is  given  by 

A1-^  =  f j  U-a-mai/is-f)'*]  =  hx  (A-22) 

Then,  the  moments  for  this  geometry  become,  for  the  absorber 
(central  region) 

&A  ji  s  Jx  ^AWX„^rt  e  A'y  (A-23) 

while  for  the  moderator  (surrounding  medium) 

+  J2  tUiSMi  I*  (&*»   +  MFC  ^fexrtj  e"*1* 


(A-24) 
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Solution  for  the  Unit  Cell  Model 
with  Isotropic  Source  Term 

In  the  unit  cell  model  it  is  assumed  that  the  z-dependence 
of  the  angular  flux  can  be  neglected,  and  that  the  symmetry  con- 
dition given  by  Eq.(A-lO)  applies.   It  is  also  assumed  that  the 
slowing  down  density  is  zero  in  the  fuel  and  constant  in  the  mod- 
erator. Thus,  the  equations  are  homogeneous  in  the  fuel  only. 

For  the  equations  in  the  absorber  and  for  the  homogeneous 
part  of  the  equations  in  the  moderator  Eq.(A-4)  and  Eq.(A-6)  will 
apply  with  X  =  0.   It  is  seen  from  Eq.(A-ll)  that  with  A  «  0  the 
R^C^cy  having  S  +  w  odd  are  completely  independent  of  the  R/M(<*.) 
with  J-+   m.  even.  However,  the  angular  flux  should  be  symmetric 
about  any  plane  perpendicular  to  the  z-  axis.   This  requires  that 

$(r,e>}^   -  fYr.r-e,  0)  (A-25) 

By  retaining  only  the  moments  with  J  +  n   even  this  condition  will 
be  satisfied. 

Applying  the  above  modifications  to  Eq.(A-ll),  and  using 
Cramer's  rule  on  the  resulting  matrix  yielded  the  following  char- 
acteristic values. 

<    *  fsili-^i-my^/s^f^tl  (A-26) 

e*  -  7 

where  g  is  given  by  Eq.(A-l8). 

Assuming  that  there  is  a  constant  isotropic  source  term  in 
the  moderator  and  no  neutron  source  in  the  absorber,  the  solution 
for  the  unit  cell  model  in  the  absorber  (central  region)  is 

i.  fr1  =  A  ^  *'M  **<*** O  ( A-27 ) 
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while  in  the  moderator  (surrounding  medium) 

iiM  =  %*»*(&fe*JeA*  +(-i?Ct&feir>J  +|fU  (a-28) 

where  the  R^,( «A)   are  given  in  Table  A-2. 

Table  A-2 

Elements  of  Solution  Vectors 
for  the  Unit  Cell  Model 
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APPENDIX  B 

Experimental  Work 

In  the  experiment  described  here,  a  steel  rod  was  inserted 
into  the  thermal  column  of  a  reactor  and  gold  "foils"  were  used 
to  measure  the  neutron  flux  in  and  near  the  rod.   This  experiment 
was  conducted  in  order  to  provide  experimental  data  for  a  test  of 
the  computer  code  described  in  Appendix  D,  to  illustrate  a  pos- 
sible procedure  for  the  determination  of  an  "equivalent  rod 
radius*,'  and  to  provide  a  rough  comparison  between  experiment  and 
theory. 

The  thermal  column  of  the  Argonaut  Reactor  at  Argonne  Nat- 
ional Laboratory,  Argonne,  Illinois  was  used  for  this  experiment. 
Details  of  the  Argonaut  Reactor,  its  thermal  column,  and  the 
"standard"  one-slab  loading  used  here  can  be  found  in  several 
references  (13,21). 

A  special  arrangement  which  was  first  used  by  L.  Seren  (17) 
is  shown  in  Fig.  B-l.   These  special  stringers  replaced  the  reg- 
ular J-10  (central)  and  J-ll  stringers  in  the  Argonaut  thermal 
column.  A  one  inch  diameter  steel  rod,  24-inches  in  length,  was 
placed  near  the  core  end  of  the  stringer  in  the  position  provided 
and  the  remainder  of  the  one-inch  channel  was  filled  with  one- 
inch  graphite  cylinders.  The  two-inch  channel  was  filled  with 
1  and  15/16  -inch  graphite  cylinders-. 

The  24-inch  steel  rod  v&s   cut  into  sections  10,  8,  and  6 
inches  long.   Slots  £  inch  deep  and  1/16  inch  wide  were  cut  dia- 
metrically into  both  ends  of  the  10-lnch  section.   The  6-inch 
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section  was  placed  at  the  core  end  of  the  stringer,  and  the  three 
sections  were  placed  end  to  end  so  as  to  give  the  appearance  of  a 
solid  rod  having  foil  slots  approximately  6  and  16  Inches  from 
the  core  end  of  the  rod  (or  approximately  13  and  23  inches  from 
the  core  end  of  the  J-10  stringer).  The  Cadmium-ratio  at  these 
two  points  is  15  and  45,  respectively  (20)v. 

Foils  and  Counting  Facilities 

Two  types  of  foils  were  used  to  obtain  the  fine  structure 
near  the  iron  rod.  One  was  a  uniform  gold  wire  with  a  32-mil 
diameter.  The  other  was  a  gold  ribbon  i-inch  wide  by  2-milB 
thick.  Both  types  of  foils  were  cut  into  strips  approximately 
5i-inches  long. 

The  counting  system  used  with  the  gold  wire  is  shown  in  Fig. 
B-2.  The  shielded  end-window  G.  M.  tube  was  a  Nuolear  Chicago 
Model  3031B,  serial  344,  with  U.  S.  Govt,  number  93541.  The 
scaler  and  power  supply  were  designated  by  U.  S.  Govt,  number 
92121.  The  central  lead  shield  had  originally  been  designed  with 
a  i  inch  diameter  opening  at  the  top  and  with  an  intersecting 
wire  channel  of  approximately  £  inch  diameter.  The  central  open- 
ing was  reduced  to  1/16  inch  by  insertion  of  a  small  lead  plug. 
The  intersecting  wire  channel  was  reduced  to  approximately  l/l6 
inch  by  means  of  plastic  tubing. 

The  gold  ribbon  was  counted  in  an  end  window  gas  flow  pro- 
portional counter.  The  counter  had  U.  S.  Govt,  number  92665,  and 
the  scaler  and  power  supply  had  U.  S.  Govt,  number  92130. 
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Experimental  Procedure 

The  5£  Inch  long  foils  were  placed  along  the  foil  grooves 
and  through  the  steel  rod  so  that  one  end  of  the  foil  was  flush 
with  the  edge  of  the  J-ll  stringer.  The  reactor  was  operated  at 
a  nominal  power  of  75  watts  for  one  hour.   The  stringers  were 
removed  with  the  aid  of  a  threaded  rod  which  screwed  into  the  end 
of  the  special  stringers.  After  removal  of  the  foils  the  iron 
rod  was  allowed  to  "cool"  overnight  before  reusing. 

The  gold  wire  was  counted  in  the  apparatus  shown  in  Fig.  B-2. 
The  wire  was  visually  positioned  with  the  aid  of  the  mounted 
scales  whose  smallest  division  was  1/16  inch.  The  portion  of  the 
wire  which  had  been  in  or  near  the  rod  was  counted  at  l/l6  inch 
intervals.  The  remainder  of  the  wire  was  counted  at  l/8  inch 
intervals.   The  wire  was  counted  for  pre- set  times  ranging  from 
2  to  8  minutes,  depending  on  the  count  rate.   The  appropriate 
counting  time  was  chosen  so  that  the  standard  deviation  of  the 
total  count  at  any  position  was  less  than  2%   of  the  total  count. 

The  gold  ribbon  was  cut  into  foils  approximately  1/20  of  an 
inch  in  length  by  taping  the  ribbon  on  graph  paper  and  cutting 
the  ribbon  at  each  division.  Every  foil  which  had  been  in  or 
near  the  rod  was  counted  and  every  third  foil  was  counted  for  the 
remainder  of  the  region.  The  foils  were  counted  in  the  gas  flow 
end  window  proportional  counter  for  pre- set  times.   The  same 
criteria  was  used  to  determine  the  counting  time  as  was  used  for 
the  gold  wire.  After  counting,  each  foil  was  weighted  on  an 
electric  balance. 
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Treatment  of  Raw  Data 

All  counts  were  corrected  for  background,  and  those  taken 
with  the  0.  M.  tube  were  corrected  for  dead  time.   The  dead  time 
correction  for  the  proportional  counter  was  insignificant.  The 
activity  was  then  corrected  for  decay  of  the  gold  nuclei  using 
the  well  known  exponential  decay  law  and  a  half  life  of  2.70  days 
(5).  Since  only  relative  activities  were  of  interest,  the  meas- 
ured activity  at  a  given  position  was  corrected  to  the  activity 
present  when  the  first  section  of  the  wire  or  ribbon  was  counted. 
The  data  taken  with  the  gold  ribbon  was  placed  on  a  milligram 
basis  by  dividing  the  activity  of  each  foil  by  the  weight  of  the 
foil.  All  counts  obtained  from  a  given  5i-inch  gold  strip  were 
then  placed  on  the  same  time  basis  so  that  relative  count  rates 
at  different  positions  could  be  compared. 

The  data  processed  in  the  above  manner  are  tabulated  in 
Tables  B-l  through  B-3.   The  deviation  listed  is  the  standard 
deviation  of  the  net  count  rate.   In  computing  the  standard  devi- 
ation for  the  activity  from  the  gold  ribbon  it  was  assumed  that 
the  standard  deviation  in  the  weights  of  the  foils  was  0.025  mg. 
This  figure  was  arrived  at  by  the  following  considerations: 
Since  the  foil  weights  could  be  read  accurately  to  0.1  mg,  the 
maximum  error  Bhould  be  0.05  mg.   If  it  is  assumed  that  the 
"maximum"  error  is  of  the  order  of  twice  the  standard  deviation, 
then  0.025  mg  is  the  standard  deviation  In  the  weight  of  a  given 
foil.  The  combined  standard  deviation  was  then  calculated  in  the 
usual  manner  (5). 
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Table  B-l 


Processed  Experimental  Data  for  1"  Rod 

Using  Gold  Wire  at  Approximately  23" 

from  Core  End  of  J-10  Stringer 


Distance  from 

Counts  per 

Distance  from 

Counts  per 

Rod  Center 

4-min 

Rod  Center 

4-mj 

.n 

(1/32) -inches 

(1/32) -inches 

-33 

12292 

t 

115 

21 

11114 

*  108 

-31 

12219 

t 

115 

23 

11755 

*  111 

-29 

12050 

* 

114 

25 

11938 

*  112 

-27 

11937 

* 

113 

27 

12039 

*■  112 

-25 

11552 

i 

111 

29 

12158 

*  113 

-23 

11159 

i 

109 

31 

12164 

*  113 

-21 

11104 

+ 

109 

33 

12327 

4  113 

-19 

10925 

£ 

108 

37 

12741 

i  115 

-17 

10504 

* 

106 

41 

12775 

*■  115 

-15 

10146 

X 

104 

45 

12808 

*  115 

-13 

9615 

L 

101 

49 

12916 

*  116 

-11 

9235 

+ 

99 

53 

12983 

*  116 

-9 

8768 

+_ 

97 

57 

13033 

*  116 

-7 

8836 

+ 

97 

59 

13406 

i  118 

-5 

8710 

^ 

96 

61 

13235 

*  117 

-3 

8438 

*j 

95 

63 

13004 

*  116 

-1 

8462 

t 

95 

65 

13190 

i  117 

1 

8435 

X 

95 

69 

13188 

*  117 

3 

8387 

X 

95 

73 

13508 

*•  118 

5 

8613 

± 

96 

77 

13546 

i  118 

7 

8783 

X 

97 

81 

13864 

*  119 

9 

8922 

X 

C7 

85 

13661 

*  118 

11 

9347 

± 

99 

89 

13257 

i  117 

13 

9744 

X 

101 

93 

13230 

t   116 

15 

10325 

X 

104 

97 

13174 

*  116 

17 

10585 

X 

106 

101 

13179 

•t  116 

19 

11046 

X 

108 

105 

13162 

*  116 
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Table  B-2 


Processed  Experimental  Data  for  lw  Rod 

Using  Gold  Wire  at  Approximately  13" 

from  Core  End  of  J-10  Stringer 


Distance  from 
Rod  Center 
(1/32) -inches 

-34 
-32 
-30 
-28 
-26 
-24 
-22 
-20 
-18 
-16 
-14 
-12 
-10 

-8 

-6 

-4 

-2 
0 
2 
4 
6 
8 

10 

12 

14 

16 

18 

20 

22 


Counts  per 
2-min 


13478 
13173 
12886 
13021 
12844 
12438 
12379 
12038 
11753 
11350 
11107 
10860 
10695 
10246 
10036 
10063 
9982 
10115 
10059 
10001 
10112 
10439 
10637 
10950 
11062 
11593 
11893 
12100 
12358 


i  118 

*  117 

*  116 

*  116 

*  115 
i  113 

*  113 
*■   112 

*  110 
i  108 

*  107 

*  106 

*  105 

*  103 

*  102 
*■   102 

*  102 

*  102 

*  102 
*■  102 

*  102 

*  104 
±  105 

*  106 

*  107 

*  109 

*  110 
±  111 

*  113 


Distance  from 
Rod  Center 
(1/32) -inches 

24 

26 

28 

30 

32 

34 

36 

38 

40 

42 

44 

46 

48 

50 

54 

58 
•  62 

66 

70 

74 

78 

82 

86 

90 

94 

98 
102 
106 


Counts  per 
2-min 


12424 
12935 
12983 
13322 
13671 
13569 
13659 
13840 
13682 
13860 
13949 
13992 
13959 
14007 
14375 
14562 

14763 
14941 

14769 
14894 
15081 
15189 
14766 
14617 
14261 
14168 
14061 
14086 


*  113 
±  115 

*  115 

*  117 

*  118 

*  118 
±  118 

*  119 

*  118 

*  119 

*  119 
*■  119 
±  119 

*  119 
±  121 

*  122 

*  122 

*  123 

*  122 

*  123 

*  124 

*  124 
±  122 

*  122 

*  120 

*  120 
±  119 

*  119 
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Table  B-3 


Processed  Experimental  Data  for  l"  Rod 

Using  Gold  Ribbon  at  Approximately  13" 

from  Core  End  of  J-10  Stringer 


Distance  from 

Counts 

3  per 

Distance  from 

Counts  per 

Rod  Center 

2-min  per  mg 

Rod  Center 

2-min  per  mi 

(1/40) -inches 

(1/40) -inches 

-50 

15391 

£ 

66 

8 

10017  *•   48 

-48 

15580 

£ 

70 

10 

10093  *  48 

-46 

15530 

A 

71 

12 

10293  *  48 

-44 

15268 

4 

65 

14 

10996  *  60 

-42 

15345 

1 

71 

16 

11552  *  49 

-40 

14982 

A 

75 

18 

11947  •*  57 

-38 

14630 

± 

62 

20 

12471  *   57 

-36 

14508 

1 

63 

22 

13269  *■  63 

-34 

14728 

± 

67 

24 

13660  *  60 

-32 

14095 

1 

64 

30 

14318  i  64 

-30 

13718 

X 

58 

36 

14796  *  67 

-28 

13375 

1 

56 

42 

15032  *  62 

-26 

13125 

± 

60 

48 

15178  ±  63 

-24 

13009 

A 

58 

54 

15374  *  70 

-22 

12734 

1 

60 

60 

15738  *  71 

-20 

12364 

1 

58 

66 

15797  *  67 

-18 

11715 

A 

51 

72 

16105  -  76 

-16 

11133 

1 

53 

78 

16307  *  67 

-14 

11056 

4 

55 

84 

16261  *  69 

-12 

10905 

± 

51 

90 

16217  *•  71 

-10 

10481 

i 

53 

96 

16024  *.  65 

-8 

10297 

1 

52 

102 

16604  *  79 

-6 

10107 

A 

52 

108 

16145  *  69 

-4 

9956 

£ 

49 

114 

16652  a  74 

-2 

9857 

A 

50 

120 

16622  *  68 

0 

9798 

X 

48 

126 

16636  *  71 

2 

9678 

A 

46 

132 

16399  ±  65 

4 

9714 

± 

48 

138 

16707  *  68 

6 

9739 

A 

50 

144 

16917  *  77 
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Moderator  Radius  and  Cross  Sections 

In  order  to  use  the  P^  approximation  to  the  Boltzmann  equa- 
tion as  described  in  the  theory,  the  absorption  and  transport, 
or  scattering  cross  sections  in  both  moderator  and  absorber,  as 
well  as  the  physical  dimensions  of  both  moderator  and  absorber, 
must  be  known.  Although  the  dimensions  of  the  absorber  were  well 
known,  the  dimensions  of  an  equivalent  free  surface  cylindrical 
moderator  to  replace  the  actual  thermal  column  and  surrounding 
shield  had  to  be  estimated. 

An  attempt  was  made  to  determine  the  equivalent  moderator 
radius  by  making  traverse  flux  measurements  in  the  Argonaut  ther- 
mal column  and  fitting  the  data  with  a  series  of  orthogonal  Bessel 
functions  of  the  first  kind  of  zero  order.  As  predicted  by  diff- 
usion theory,  the  equivalent  radius  would  then  be  the  point  where 
the  resulting  curve  goes  to  zero. 

The  use  of  Just  the  first  term  of  the  series  yielded  an 
equivalent  radius  of  approximately  30  inches,  which  is  the  actual 
dimension  in  the  direction  of  the  measurements  and  also  in  the 
direction  along  the  length  of  the  5£  inch  fine  structure  foils. 
The  data  was-  so  erratic  however,  that  higher  harmonics  yielded 
unbelievable  equivalent  radii. 

Inasmuch  as  the  flux  shape  in  and  near  the  absorber  should 
be  a  slowly  varying  function  of  moderator  radius,  for  large  radii, 
and  since  the  first  harmonic  yielded  a  radius  close  to  the  actual 
radius,  the  equivalent  moderator  radius  was  taken  as  30.0  inches. 

No  direct  information  was  available  on  the  thermal  column 
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cross  sections,  however  a  measurement  was  available  which  shows 
the  exponential  decrease  of  the  flux  along  the  length  of  the 
thermal  column,  so  the  absorption  cros-s  section  was  calculated 
from  the  equation  obtained  from  the  P]_  approximation 

where  y  s  -  tyti   and  z  is  the  direction  along  the  length  of  the 
thermal  column.  The  scattering  cross  section  was  taken  to  be 
0.285  cm"1  at  a  graphite  density  of  1.6  g/cm^  (16),  y  was  deter- 
mined by  Springer  ('20)'  to  be  0.0370  cm""1,  and  a  and  b  were  taken 
to  be  the  actual  physical  dimensions  of  24  and  30  inches.  The 
absorption  cross  section  of  the  moderator  was  thus  calculated  to 
be  2.53x10"^  cm"1. 

The  scattering  cross  section  of  the  iron  rod  was  taken  as 
0.947  cm"1  (16).  At  the  time  of  this  analysis  the  iron  absorp- 
tion cross  section  was  not  known,  and  the  value  of  0.259  cm"1 
waa  assumed  (20).   It  was  later  found  that  the  actual  absorption 
cross  section,  as  determined  by  several  measurements  at  Oak  Ridge 
National  Laboratory,  is  0.229  cm""1. 

The  total  flu::  should  not  be  a  rapidly  varying  function  of 
the  rod  absorption  cross  section  and  a  similar  comparison  between 
theory  and  experiment  would  be  expected  if  the  correct  cross 
section  had  been  used. 
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APPENDIX  C 

Description  and  Explanation  of  IBM-1620 
Computer  Program  used  to  Calculate  the 
Angular  Flux  at  the  Interface 

This  program  calculates  the  angular  flux  distribution-  at 
the  interface  between  two  media  by  using  the  P^  approximation 
to  the  Boltzmann  equation  for  a  system  having  cylindrical  geo- 
metry and  exponential  flux  dependence  in  the  z-directlon.  The 
program  was  written  in  FORTRAN. 

The  angular  flux  for  this  approximation  is 

f(V,  },e>t(p)    =XtL(r,y>(lJe,*)  (0-1) 

where  fprK  *  (-l)mfjmt   and  the  fpM  are  as  follows: 
In  the  fuel  (central  region), 

&*  ft  P  *  Ja  &/<>  IJ«+tr)  £**  ( c-2 ) 

In  the  moderator  (surrounding  region), 

+  ik/fctffcl-teirt  *  (-OX  /C^zrxK*2*       (°"3) 
The  Rj>„,  *  ,  and  £  are  given  in  Appendix  A.  The  A.,  Blf.  and 
CjL  can  be  determined  from  two  sets  of  boundary  conditions  with 
this  program.  One  set,  which  is  called  Case  2,  matches  all  mom- 
ents except  the  f^Q  and  f^2  at  the  interface,  and  uses  the  equa- 
tions obtained  by  using  the  P10,  Pni  P31,  and  P33  in  Marshal "s 
boundary  conditions  at  the  moderator  free  surface.  The  other  set 
called  Case  4,  matches  all  moments  except  the  f2Q  and  f,0  at  the 
interface,  and  uses  the  equations  obtained  by  using  the  p.-,  P,.. , 
P"32»  and  P-53  in  Marshak's  boundary  conditions. 
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The  unknowns  A^,  Bj_,  C*,  and  A.  are  determined  with  the  use 
of  the  matrix  equation  TX  «  0,  where  X  is  a  column  matrix  con- 
taining the  A£,  B*,  and  Cif  and  T  is  a  square  matrix  containing 
the  coefficients  as  given  by.  the  boundary  conditions.  The  cor- 
rect value  of  A.  is  taken  to  be  the  one  which  makes  the  magnitude 
of  the  determinant  of  T  the  smallest.  A.,  is  arbitrarly  set  equal 
to  1.0  and  the  other  x^  are  determined  from  the  first  11  equa- 
tions. The  determinant  of  T  and  the  x,  are  calculated  with  the 
aid  of  a  special  subroutine,  CRAM(X),  which  uses  18  digit  arith- 
metic and  performs  a  Crout  reduction  with  an  auxiliary  matrix 
using  the  method  described  by  Hildebrand  (10 ).  The  value  of  the 
determinant  is  calculated  by  multiplying  together  the  diagonal 
elements  of  the  auxiliary  matrix. 

Special  subroutines  are  used  to  calculate  the  values  of  the 
zero  and  first  order  Bessel  functions  and  the  recursion  relatione 
are  used  to  calculate  all  higher  order  functions  except  the 
J2(£i2r)  and  Jj(&i.zr)   evaluated  at  the  interface.  The  argument 
is  small  in  this  case  and  the  first  two  terms  of  the  series  expan- 
sion are-  used. 

The  Rj>w  are  stored  in  the  C0EF(I,j)  matrix,  the  <*<;,  and  £,• 
are  stored  in  TERM(J),  the  &±   and  bA  of  Table  A-l  are  stored  in 
ROOT(J),  and  the  N^'a^  and  NjCbi)  of  Table  A-l  are  stored  in 
OFR('J). 

The  proper  sequence  for  loading  the  input  data  is  shown  in 
the  source  program.  The  meaning  of  the  symbols  is  given  in 
Table  C-l. 
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Table  C-l 


Input  Data 


Symbol 


Explanation 


MM38 
MM12 


mm 

GAM0(1 

GAM0(2 

GAMl(l) 

GAMl(2) 

SIGF 

RF 

SIGM 

RM 

CRIT1 

DL 

AL 


■  0  (This  was  used  In  program  testing) 

=  0  When  A  unknown 

s  1  When  A  is  known  from  a  previous  calculation 

and  only  angular  flux  is  deBired 
=  1  For  Case  4  "boundary  conditions 
»  2  For  Case  2  boundary  conditions 
la/l  in  fuel 
Ia/x  in  moderator 
Tir/l   in  fuel 
Ztr/l  in  moderator 
I  in  fuel  (cm-1) 
Rod  radius  (cm) 
I  in  moderator  ( cm~l ) 
Moderator  radius  (cm) 
Convergence  criteria  for  X      (1.0x10"^) 
Increment  in  A 
Trial  \ 


Most  of  the  output  is  self  explanatory,  however  all  input 
data  is  printed  out  immediately  after  being  read  in,  and  this 
output  is  unlabled.   The  dimensions  on  all  output  data  is  in  CGS 
units.  The  meaning  of  all  symbols  and  words  used  in  the  output 
is  given  in  Table  C-2. 

Approximately  one  and  one-half  minutes  is  required  for  one 
calculation  of  SUM.   If  the  original  estimate  of  A  is  reasonably 
close  to  the  correct  value,  the  interpolation  converges  after 
about  three  trials  (approximately  5  minutes).  A  rough  first 
estimate  of  A  can  be  obtained  from  the  ^   approximation  for  a  one 
region  medium. 
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Table  C-2 
Output 


Symbol 


Explanation 


DELTA  LAMBDA 
R(I+1)- 

PHI 

P 

THETA 

SUM 

SLOPE 

LAMBDA 

IN  FUEL 

IN  MODERATOR 

FLUX  INYfARD 

ANGULAR  FLUX 


«  0    (Difference  of  new  A.  and  old  \  Is  zero, 
program  continues) 
Values  are  listed  in  the  following  order: 
A2,A3,A4,C1,C2,C3,C4,B1,B2,B3,B4 

(note  that  A*-^  =  1.0) 
0   (radians) 
Angular  flux 
&  (radians) 

Value  of  determinant  of  T 
Change  in  SUM  divided  by  change  in  A 

A 

Angular  flux  evaluated  using  f,„   of  Eq.(C-2) 
Angular  flux  evaluated  using  f tm   of  Eq.(C-3) 

AVERAGED  OVER  THETA    J  ffr,},e,f>)<S"i0<J0 


The  sense  switches  do  not  alter  the  program  when  in  the  off 
position.  The  changes  which  occur  when  they  are  in  the  on  posi- 
tion are  given  in  Table  C-3« 

Table  C-3 
Sense  Switches 


Switch 


Operation  When  Switch  On 


1 
2 

3 
4 


Not  used 

Causes  SUM,  SLOPE,  and  LAMBDA  to  be  typed  after 

each  trial 

Causes  CRAM(X)  subroutine  to  type  out  T  matrix 

Not  used 
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LOGIC     DIAGRAM 
FOR 
APPENDIX      C 


Start: 

Read  in 
Program 


Read  in 
Data 


Calculate  ROOT(J),  CFR(J) 
Set  MMM  =  III  =JJKK  =  0 
Set  ALI  =  2xAL 


Form 
TEST  =  ALI-AL 


Calculate 
BESEL(I, 
and  COEFU, 


I 


,j,  -(^wls  :est/a vr 


Set 
TEST=-TEST 


No 


Form   T    Matrix 
(Depends  on  MM) 


Calculate  and 
Print    X(J) 
Set  JJKK=0 


Yes 


(Is   JJKK=Q) 


Is    MMI2=0 


No 


Form 

Augumented 
T    Matrix 


Print 

SUM.AL 

SLOPE 


Print 

FORMAT 

601 


<i 


IF  SENSE 
SWITCH  2 


Calculate 


SUM  = 


Calculate 
fIO'f2l  »fIO  »flh 
f22'f33,f30>SSR 


Tl 


ON 


Print    SUM 
SLOPE,  AL 


Calculate 
SLOPE 


OFF 


Set   SUM  I  =  SUM, 

ALI=AL 

AL  =  AL-SUM/SLOPE 


(Is    1 1  J  =0 


Print    SUM,  Set    111=  I, 
ALI=AL,   SUMI  =  SUM 
AL=AL+DL 


Calculate 
f3l'f20 


(If    MM- 


Set    JJKK=0, 


f30=  SSR 


(  Is  JJKK=C? 


(Is    JJKK  =  d 


Calculate    f- 
in  moderator 


Calculate 
in   fuel 


32 


Calculate    f 
in    fuel 


20 


Calculate    f20 
in    moderator 


Calculate    and 
Print   Inward 
Flux 


Calculate  and 
jr      Print 

/f  sin0d0 
o 

I 


Calculate   and 
Print  f(e|t4>.) 
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iterate:  for 

AT  BOUNDARY* 


LAMBDA  AND  CALCULATES  ANGULAR  DISTRIBUTION 
BOTH  IN  MODERATOR  AND  IN  FUEL 


9 
10 

11 

30 
90 

91 

93 

9  2 

601 


110 

112 
111 


DIMENSION  T(12»13)»X(12) .ROOT (8) .TERM (8 ) »ARG<20) 

DIMENSION  BESEL(20.4)  »GAMO(2)  .GAMK2)  »COEF(  10  ,  12  )  >  CFR  (  6  ) 

FORMAT ( E 16 . 7 » F 1 6 • 7 ♦ E 16 . 7 • E 16 . 7 ) 

FORMAT (F 16. 7/) 

FORMAT  (f-'l  6.  7,  El  6.  7,  El  6.  7/) 

FORMAT! I5»I5»I5»I5tI5»I5) 

F0RMAT(E16.7»E16.7,E16.7,F16.7/ ) 

FCRMAT(E16.7»E16.7/) 

F0RMAT(E14.7,I6.E14.7) 

READ  '+.MM38.MM12.MM 

4.MM33.MM12.MM 

1»GAM0{  1)  ,GAM0(2)  ,GAMK  1  )  .GAMK2) 

1.GAMOI  1)  ,GAMO(2) .GAMK 1) .GAMH2) 

1»SIGF»RF»SIGM.RM 

l.SIGF,RF»SIGM»RM 

1.CRIT1 .DL.AL 

3.CRIT1 .DL.AL 


TYPE 

READ 

TYPE 

READ 

TYPE 

READ 

TYPE 

J  =  l 

K  =  l 

K1=J+1 

A=l.  +GAMO(K)*(.8  +  27.*GAMl(K)/35.) 

ROOT! J)=  3  5.*A*( 1,-SORT(1.-108.*GAMO(K)*GAM1(K)/(35.*A**2)  )  J/18. 

ROOT!  J+l)=- 35.*A*(l.+SORT(  1.-108.  *GAMO(K)*GAMKK  )/(  3  5.*A**2)  )  )/18. 

ROOT( J+2) =35.*GAM1 (K )/ !8.*GAM1 (K)+7. ) 

ROOT! J+3) =7. 

DO  9  I=J,K1 

CFR( I )  =  ( ROOTf  I )-GAM0(K)*< 3 .*GAM 1 ( K ) + 7. )  )/(ROOT( I )*(ROGT(  I )  -7. )  ) 

IF( J-l )3»10,11 

J  =  5 

K  =  2 

GO  TO  12 

MMM  =  0 

JJKK  0 

AL1=2.*AL 

II  1=0 

KKK-0 

TEST=AL1-AL 

IFITEST )91 .92,93 

TEST=-TEST 

IFITEST/AL 

PRINT  601 

FORMAT! 17H 

GO  TO  316 

5  L  =  ,*,  l  *  •  *  ? 

no  i  i i  j=\ ,n 

I  FU-5)  112.111  .112 
TERM! J)=SQRT(ROOT(J)-SL) 

CONTINUE 

TERM! 5 ) ^SQRT ( SL-ROOT ( 5 ) ) 

DC    113    J=l  ,': 


CRITl)316»9'n94 
DELTA    LAMBDA    - 


0/) 
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202 
11? 

12  4 
118 


119 
121 


116 


125 


120 

122 

117 
123 

114 

206 


13! 


130 


132 


ARC  (J) 
ARG ( J+ 

A  R  G  (  J  + 

*\  n  r  ;  j  .l 

ARG  (J-1 
CO    11 A 
IF< J-5 
BESEK 
BESEK 
A=    -1. 
B=    -A 
GO    TO 
BESEK 
BESEK 
GO    TO 
IFU-9 
BESEK 
BESEK 
A  =  l. 
B^A 
GO    TO 
BFSFL! 
BESELI 
IF( J-9 
BESEK 
BESEK 
GO    TO 
A  =  l. 

GO    TO 
IF(  J-l 

ifu-i 

BESEK 

BESEK 

CONTIN 

A  =  -l. 

B-A 

C---A 

L  =  l 

M=2 

K  =  L 

DO    130 

COEF( 1 

COEF(2 

COEFI3 

C0EFI4 

C0EFI5 

CCEF(6 

C0EF(7 

COEFd 

COEF(8 

GOEF(9 

IF(M-2 

L  =  5 


=  T  E  R  M  (  J  )  *  S I G  F  *  R  F 

'.  )-TERM(  J+4)*S  IGM*RF 

C  s--h?c-  (  JA-I  ) 

12 ) - TF  P v ( J+ 4 ) *  S I GM*RM 
!M=APG(J  +  12) 

J-l  *20 
)  1 1 5  » 1 1 3  » 1  \  o 
J*1)=TZER( »RG( J) ) 
J»2)  =  IONE( ARG( J)  ) 


123 

J*1)=YZER(ARG( 

J»2)=YONE(ARG( 

120 

) 121 » 116*122 

J.1)=KZFR(ARG( 

J»2)=KONE(ARG( 


123 

J»1)=JZER(ARG( 

J*2)^JONE(ARG( 

)12O,125»120 

J»3)=( ARG (J >** 

J,4)=< ARG (J)** 

114 

124 

3  )  U  5  ,  11  fi  ,  l  1  7 

7)121tll6tll5 

J»3)=2.*A*BESE 

J,4)=4.*A*BESE 
UE 


J)  ) 
J)  ) 


J)  ) 
J)) 


J)  ) 

J)  ) 

7  )*( l.-ARG( J)**2/12. )/8. 
3)*(1.-APG( J)**2/16.)/48. 


LIJ.2) /ARG( J) 
KJ,3)/ARG(  J) 


+  B*BESEKJ»1) 
+  B*BESEKJ,2) 


J  =  L»M 
»J)  =A 

»J) =B*6.*TERM( 
, J) =A*3.*GAM0( 
.J) =A*(3.*SL-R 
» J) =C*TERM( J)* 
*J) =-B*3.*(SL- 
. J) =3.*AL*COEF 
C* J ) =A*3.*AL*( 
.J) =C*TERV( J)* 
*  J) =B*9.*TFRM( 
)  132,  K?,  131 


i)*AL*CFR(J) 
K)*AL/ROOK  J) 
OOT( J) )*CFR(J) 
COEF(3.J)/AL 
ROOK  J  J  )*CFR(  J) 

(  6  ,  J  ) 

5.*SL-3.*ROOK J) )*CFR(J)/5. 

CCEF(6»J)/2. 

J)*(5.*SL-ROOT(J) )*CFR(J)/1C. 


80 


136 


131 
135 


13 '4 


142 


K=2 

A  =  -A 

C  =  -C 

M=L 

GO    TO    133 

IFIM-M  135,135»134 

L  =  6 

3^-B 


GO    TC 
K  =  l 
A  =  -l. 
J  -  3 

conn  i 

C0EF(2 

cocr(3 

CCEF(4 
CCEFI5 
C0EF(6 
C0FF!7 
C0EF(8 
C0FF(9 
COEFt 1 
CCEFd 
CCEF(2 
C0EF(3 
C0EF(4 
C0EF(5 
C0EF(6 
CCEF17 
C0EF(8 
C0EF19 
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J)=0. 

J) =-A*5.*GAMl IK)*(2.*SL-R00T( J) )/( TERM( J)*RCOT ( J) ) 

I)    .A 


J)  =A*5.*GAM1  IK.)  *AL/POOT  (  J) 
J) =-A*AL/TERM(J) 
J)=-C0EF(4iJ) 

J) =-A*5.*GAMl (K)*(3.*SL-R00T( J) )/ROOT( J) 
J)=-TFRMU)*C0EF(4»J)/2. 

J)=-A*GAMi  (K)*AL*(  15.*SL-11.*ROOT( J)  )  /  (  2  .*ROOT  (  J  ) J' TERM  (  J  }  ) 
»J >=A*GAM1 (<)*( 5.*SL-R00T( J) )/ROOT( J) 
J+l)  = 0.0 

J+l) =-A*2.*AL/TERM( J+l ) 
J+D-0.0 
J+1)=A 
J+l)=C.O 

j+1 )=-A*(SL+ROOT( J+l ) )/(SL-Rn0T( J+l) ) 
j+1 )=-a*al*( ^.*£L-R0nT( J+L ) ) / ( SL-ROOT ( J+l ) ) 
J+l )^A*(SL+ROOT( J+l ) ) /( 2.* TERM (j+l j ) 
J+l  )--A»(3.*SL-R00T( J  +  l )  )/(2.*TERM( J  +  l)  ) 
CCEF(10*J+1)-A*AL 
IF(  J-  3)  141,140,141 
14C    J  =  7 
A- -A 
K^2 

GO    TO    142 
14  1    00    146    K=5*8 

00    147    J=I » 1C 
147    C0EF(J»K+4)=COEF(J»KJ 

I F(K— 5 ] 146,146*148 
14e    CO    149    J=2,8,3 
149    COEF( J,K) ~-COEF( J  >  K  ) 
C0EF(9,K>=-C0FF(9»K) 
146    CONTINUE 
211    DO    150    J=l  ,12 

DC    151    KM, 3.2 
151    T(K»J)=C0EF(K»J)*BESEL(J»1 ) 
IF(MM-1 )713. 714*71° 

714  T  (4  ,J)=C0EF(9.J)*BFSEL( J.2) 
GO    TO    715 

7  13    T('4,J)=C0EF  (4,J)*BESEL<  J.l  ) 

715  DO    15?    K.=  2.5»3 
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152  T(K  »J)=CO'r'r(<»J)*"f:SEL  (  J»2  ) 

T  (6,J)=COE<r(6*J)*BESEL(J»3) 

I  F(MM-1  !  73C  ,  730»731 
73  0  T(7,J)=C9EF( 7,J)*BESEL(J»3) 

GO  TC  150 
731  T(7»J)=C0EF(9»J)*BE5ELU»2) 
150  T(8,J)=CCEF(e*J)*8E5EL(J»4) 

DO  155  J=9»12 

DC  155  K=1.4 

155  T(J,K)=0.0 

DO  156  J=5»12 

T(9»J)  =C0EF(1»J)*SESEL( J+8.1)  +  2 . *COEF ( 5 , J ) *BESEL ( J+8 , 2 ) /3 . 
T (9»J)  =  T(9,J)-(C0EF(4.J)*8ESEL( J+8 , 1 ) -COEF ( 6 , J ) *BESEL ( J+8, 3) )/8. 
IF(MM-1 }720»721,720 
72  0  T( 10,J)=?.*COEF(^  » J ) *BESEL ( J+8 » 1 ) /3.  +COEF ( 2 » J } *BESEL ( J  +  9  ♦  2 ) / 4 . 
GO  TO  722 

721  T(10,J)=-C0EF(6,J)*BESEL( J+8,3)/48.  +8 . *COEF ( 9 , J ) *BESEL ( J+8 ,2 ) /2 1 . 
T(10»J)=T(10»J)  +(COEF(  1,J)M.  +7  .  *COEF  (  4  ,  J  )  /  16  .  )  *BESEL  (  J  +  8  , 1  ) 

722  I F(MM-l) 723. 723*724 

72  3  T(  ll»J)=COEF(2»J)*BESEL(J+8.2)/8.  +4.  *COEF  (  7  »  J  )  *6ESEL  (  J  +  8  .  3  J  /35  . 

GO  TO  725 
72 4  T(  ll»J)=-COEF(6.J)*3ESEl(  J+8»3)M8.  +8  .  *COEF  (  9  » J  )  *BESEL  (  J  +  8  »2  ) /21. 

T(ll»Ji=T(ll*J)  +(C0EF(l,J)/4.  +7.*COEF(4,J )  / 16. )*BESEL ( J  +  8 » 1 ) 
72  5  T(12»J)=-COEF(l  ,  J  )  *BESEL  (  J^-8  ,  1  )  /4  .  +8 . *COEF ( 8 , J ) *BESEL ( J+8 .4 ) /35. 

T(12tJ)=T(12tJ)  +COEF(4»J)*BESEL( J+8»l)/16. 

156  TU2.J)=T(12*J)  +3.*COEF(6.J)*BESEL(J+8.3)/16. 

160  IF( JJKKJ1 160*1160 »944 
1160  IF(MM12)161»161»3C'2 

161  IF(MM38)902»PU2.90l 
902  SUM=CRAM( -12. ) 

GO  TO  299 
901  X( 1  )=CRAM(  12. ) 

PRINT  610 
610  FCRMAT15X7H  R(I+1)/) 

DO  935  J=l ,10 
93^  PRINT  1,X(J) 

PRINT  2»X(  11  ) 

JJKK=1 

GO  TO  211 
944  DC  945  J=l ,8 
94  5  ARG(J)=  -T(J.l) 

DO  946  J=l,8 

DO  94  6  K=l,3 
946  ARG(J)=ARG(J)  -X ( K ) *T ( J  ,K+ 1  ) 

F00=ARG( 1 ) 

F21=ARG<2) 

F10=ARG(3) 

F11=ARG(5  ) 

F22=ARG(6  ) 

F33=ARG{ 8) 

SSR=-COEF( lO»l)*BFSEL(ltl) 

F30=C. 

DO  1618  J =1,3 
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1618  S5R=SSR  -X(J)*C0EF(10.J+] )*BESEL< J+l»l) 
DO  1619  J=4»ll 

1619  F3O=F30  +X ( J ) *COEF ( lu , j+l ) *BESEL ( J+l » 1 ) 
1513  IFIMM-1 )8, 500, 501 

500  F31=ARG(4) 
F32=ARG(7) 
IF(JJKK)1500, 1500,1501 

1501  F20=0. 

DO  1502  J=4,ll 

1502  F20  =  F20  +X ( J ) *COEF ( 4 , J  +  l ) *BESEL ( J  +  l , 1  ) 
GO  TO  502 

1500  F20=  -C0EF(4,1)*BESEL(1,1) 
DO  1503  J=l»3 

1503  F20=F20  -X ( J ) *COEF ( 4 , J+l ) *BESEL ( J+l , 1 ) 
GO  TO  502 

531  F31=ARG(7) 
F2C=ARG(4) 
IF{ JJKK)1521 ,1521*1522 

1522  F32=0. 

DO  1523  J=4,ll 

1523  F32  =  F3?  +X ( J ) *COFF ( 7 , J+l ) *BESEL ( J  +  l ,  3  ) 
GO  TO  502 

1521  F32  =  -COEF(7»l)»BESEL(l,3) 
DO  1524  J  =  l,3 

1524  F32=F32  -X ( J ) *COEF ( 7 , J+l ) *BESEL ( J+l , 3 ) 

502  IF(JJKK)503, 503,504 

504  PRINT  698 

698  FORMATC20X13H  IN  MODERATOR/) 
GO  TO  505 

503  PRINT  699 

699  FORMAK20X8H  IN  FUEL/) 

505  PRINT  622 

622  FORMAT(7X4H  PHI.13X2H  F.12X6H  THETA/) 

THETA=0. 

F=F00  +F10  +F20  +F30 

PRINT  666, F, THETA 
666  FORMATI3X13H  ALL  VALUES   ,  E16.  7  ,  E  16.  8  ) 

DO  350  1=1 ,5 

PHI=0.0 

THETA=THETA+, 5235988 

CY=COS( THETA) 

SY=SIN< THETA) 

IF( 1-3)351 ,532,351 
5?2  CY=0. 
351  DO  ^50  J=l ,1 1 

CPHI=COS(PHI ) 

CPHI2=COS(2.*PHI ) 

f PHI3=COS(3.*PHF ) 

IF(  J-6J37.1  ,372*371 
372  CPHI=0.U 

CPHI3=0.0 
371  F=F00  +F10*CY  -F11»SY*CPH1  +  F20* ( 3 . * ( C Y**2 ) -1 . ) / 2 . 

F  =  F-F21*SY*CY*CPHI  +  F22* ( SY**2 > *CPH I  2/2 . 
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F=F  +F30*(5.*(CY**3)-3.*CY)/2<  +F32*( SY**2 )*CY*CPHI 2/2, 

F=F-F31*SY*(5.*(CY**2)-1. )*CPHI/2.  -F33* ( SY**3 ) *CPHI 3/2. 

IF(J-1)359,359,358 
359  PRINT  1 »PHI »F,THETA 

GO  TO  35U 
358  PRINT  1 »PHI  tF 
35  0  PHI = PHI  +.3141593 

THETA=3. 1415927 

f-=  FOO  -F1U  +  F20  -F3G 

PRINT  666,h,THFTA 

PRINT  531 
631  FORMAT (YX4H  PHI,3X33H  ANGULAR  FLUX  AVERAGED  OVER  THFTA/) 

PI  =  3.1415927 

PHI  =  0.0 

DO  766  I  =1,11 

CPHI=COS(PHI ) 

CPHI2=COS(2.*PHI ) 

CPHI3=COS(3.*HHI ) 

I F ( I  -6 )  382,381,382 

381  CPHI=C.O 
CPHI3=0.0 

382  FA=2.*F00  -P I*F11*CPHI /2.  +2 .*F22*CPHI 2/3 . 
FA=FA-PI*F31*CPHI/16.  -P I *3 . *CPHI 3*F33/ 16 • 
PRINT  1, PHI, FA 

766  PHI  =  PHI  +.3141593 

FIN=PI*(2.*F00  +F11  +F31/8.  -F33/8.) 

PRINT  638»FIN 
638  F0RMAT(E16.7.14H  =  FLUX  INWARD/) 

IF(JJKK)8,8,1512 
1512  JJKK=0 
F30=SSR 
GO  TO  1513 

299  IF( I  I  I ) 300,300,304 

300  111=1 

PRINT  602,  SUM 
602  F0RMAT(E16.7,6H  =  SUM/) 

AL1=AL 

SUM1=SUM 

AL=AL+DL 

GO  TO  90 
304  SLOPE=( SUM-5UM1)/(AL-AL1) 

IF(SENSE  SWITCH  2)330,313 
330  IF(KKK)314»314,315 
31  '\    KKK.=  1 

PRINT  6'V, 

315  PRINT  3, SUM, SLOPE, AL 
313  SUM1=SUM 

AL1=AL 

AL=AL-5UM/SLOPE 
GO  TO  90 

316  PRINT  604 

604  FORMAT(8X4H  SUM,9X6H  SLOPE, 11X7H  LAMBDA/) 
PRINT  3, SUM, SLOPE, AL1 
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352 
317 


319 


319 


DO  317  J=] ,12 

T( J,13)=  -T( J,l ) 

DO  318  J=l.ll 

DO  318  1  =  1  ,1! 

T( J.I )=T( J,I+1 ) 

DO  319  J=l»ll 

T(Jtl2)=0. 

T(12»J)=0. 

MM38=1 

CC  TO  161 

END 
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APPENDIX  D 

Description  and  Explanation  of  IBM-1620 

Computer  Program  used  to  Determine  an 

Equivalent  Radius 

This  program  uses  an  interpolation  procedure  to  calculate 
the  fuel  rod  radius  which  makes  the  sum  of  the  weighted  squares 
of  the  residuals  "between  experimental  and  theoretical  values  of 
the  total  neutron  flux  a  minimum.  The  program  was  written  in 
FORTRAN. 

The  total  neutron  flux  in  a  two  region  cylindrical  medium 
is  described  by  the  P^  approximation  to  the  Boltzmann  equation 
for  neutron  transport.   The  total  flux  is  proportional  to  fQQ 
and  the  general  fiM   can  be  written  in  the  form 

fi-Ay*  =  j^-Oi^  (D-l) 

where  in  the  fuel  (central  region) 

G,>p  '-  iA^tJ^IJ^Xr)eKiy  (D-2) 

and  in  the  moderator  (surrounding  region) 

f  i^erf^1^**  +  ufctu$idex*2*         (D"3) 

The  R^,  vn,  and  <s„  are  given  in  Appendix  A.  The  aJ,  B1, 

"■      n* 

and  Cn,  in  terms  of  aJ  s  1.0,  and  the  A*  are  calculated  numer- 
ically in  this  program  by  the  method  described  in  Appendix  C  and 
in  the  Discussion.  For  boundary  conditions,  all  f  Jm  except  f 
and  f32  are  matched  at  the  interface,  and  the  equations  obtained 
by  using  P1Q,  Pn,  P31,  and  P33  in  Marshall's  boundary  conditions 
at  the  moderator  free  surface  are  used. 
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By  considering  only  the  radial  variation  of  the  flux  at  the 
points  r,,  and  defining 

the  sum  of  the  weighted  squares  of  the  residuals  can  be  written 


as 


E  =  i^E|D,^-$.f  (D-5) 

where  Wi  is  the  weighting  function  and  $y  is  the  experimental 
value  of  the  total  flux  at  the  point  r*. 

In  this  program  an  initial  estimate  of  the  rod  radius  is 
used  to  calculate  the  F^*,  the  D^  are  determined  by  the  least 
BquareB  procedure,  E  is- calculated,  the  rod  radius  is  changed, 
and  the  process  is  repeated.  After  three  values  of  E  have  been 
calculated,  new  values  of  the  rod  radius  are  predicated  by  a 
second  order  polynomial  of  E  as  a  function  of  rod  radius.  This 
process  continues  until  the  change  in  the  rod  radius,  divided  by 
the  rod  radius,  is  less-  than  a  specified  precision. 

The  program  Is  arranged  so  that  it  may  be  taken  off  the 
computer  after  any  given  A  has  been  determined  or  after  E  and  a 
new  rod  radius  have  been  calculated,  without  destroying  the 
mechanism  which  has  been  set  up  to  predict  new  values  of  A.  and 
rod  radius.  The  previous  values  must  be  used  as  input  in  the 
succeeding  calculation  however. 

The  proper  sequence  for  loading  the  input  data  is  shown  in 
the  source  program.  The  meaning  of  the  symbols  is  given  in 
Table  D-l. 


87 


Table  D-l 


Input  Data 


Symbol 


Explanation 


NMOD 

MMM 


GAM0(1' 

GAM0(2 

GAMlll' 

GAM1(2' 

SIGP 

SIGM 

RM 

CRIT1 

CRIT2 

FACTR 

DIV1 

DLTRF 
EERR(l) 

EERR(2)V 

RRF(1' 

RRF('2| 

phi(j; 

K 
W(J) 

A 

NC 
RF 
ADD1 
DL 

AAL(J) 
AALO(J) 
SLOPE( J 
INDEX( J 


Number  of  data  points 

=  0  when  no  previous  value  of  E  is  known 

when  one  previous  value  of  E  is  known 

when  two  previous  values  of  E  are  known 

in  fuel 

in  moderator 
in  fuel 


-  1 
=  2 

ltr/Z 


Ztr/z    in  moderator 

t  in  fuel   ( cm-1 ) 

t  in  moderator  (cm~l 

Moderator  radius  (cm, 

Convergence  criteria  for  A  (l.OxlO-8) 

Convergence  criteria  for  rod  radius 

Factor  used  to  convert  radial  position  of  data 

points  into  CGS  units 

Factor  used  to  convert  radial  position  of  data 

points  into  CGS  units 

Increment  in  rod  radius 

Next  to  last  value  of  E  if  MMM  -   2 

Last  value  of  E  when  MMM  =  1 

Last  value  of  E  when  MMM  =  2 

Rod  radius  corresponding  to  EERR(l) 

Rod  radius  corresponding  to  EERR(2j 

Data  points  in  order  from  center  out 

Any  integer 

Reciprocal  of  square  root  of  weighting  factor, 

in  order  from  center  out 

Radial  position  of  data  points  in  order  of 

increasing  r  (any  units) 

Number  of  harmonics 

Trial  rod  radius  (cm) 

1.0x10-7   (used  to  increment  X  slightly) 

Increment  in  X 

Trial  Aj  in  order  of  increasing  J  up  to  j  n  NO 

cTJ'n£w^P^evioUS  calculation,    in  order  of  Increasing   1 
SLOPE(J)    from  preceding  calculation 
*  0     if  Aj   unknown 
»  1     if  Aj   known 
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All  of  the  output  is  labled,  however  some  of  the  Input  data 
is  printed  out  immediately  after  being  read  in,  and  this  output 
is  unlabled.  The  dimensions  on  all  labled  output  is  in  CGS 
units.  The  meaning  of  all  words  and  symbols  used  in  the  output 
is  given  in  Table  D-2. 

Table  D-2 
Output 


Symbol 


Explanation 


DELTA  LAMBDA  =  0 
DELTA  RADIUS  -  0 


BUM 
SLOPE 
LAMBDA" 
R(ltl)' 


COEFICIENTS 
CALCULATED  FLUX 
RADIUS 

ROD  RADIUS 
ERROR 


(Difference  of  new  A  and  old  A  is  zero, 

program  continues) 
(Difference  of  new  rod  radius  and  old  rod 

radius  is  zero,  problem  is  finished  and 

program  goes  back  to  start) 
Value  of  the  determinant  of  T 
Change  in  SUM  divided  by  change  in  A 

Values  are  listed  in  the  following  order: 
A"2,A3,A"4,  cl»  c2»  c3»  C4,B1,,b2,B3>  B4 
(note  that  An  =  1.0) 

Di  of  Eq.(D-lT 

Total  neutron  flux  at  r. 

Radial  position  of  CALCULATED  FLUX  and  of 

experimental  data  points 

Fuel  rod  radius 

E  (Least  squares  error) 


If  the  initial  estimate  of  A  is  reasonably  close  to  the 
correot  value,  approximately  5  minutes  of  computer  time  is 
required  to  determine  the  correct  \   and  the  corresponding  x,  for 
each  harmonic  used.  For  37  data  points,  approximately  5  addi- 
tional minutes  is  required  to  calculate  the  error  and  a  new  rod 
radius  after  /  has  been  determined.  An  additional  5  minutes  is 
required  for  each  additional  harmonic. 
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The  sense  switches  do  not  alter  the  program  when  In  the 
off  position.  The  changes  which  occur  when  they  are  In  the  on 
position  are  given  in  Table  D-3. 

Tahle  D-3 
Sense  Switches 


Switch  Operation  When  Switch  On 


1  Punch  R(T-KLV 

Punch  RP,  CALCULATED  FLUX,  and  RADIUS 

2  Print  SUM,  SLOPE,  and  LAMBDA'  after  each  trial 

3  Causes  CRAM(X)  subroutine  to  type  out  T  matrix 

4  Interpolation  for  rod  radius  is  bypassed, 
new  values-  of  RF  must  be  read  in 
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Start: 

Read  in 
Program 


LOGIC     DIAGRAM 

FOR 

APPENDIX    D 

MS. 


m. 


Begin  Data 
Read  in 


Read   RF 

Calculate     ROOT(J),CFR(J) 
Set   INDXI(J)=0,  RFI=2XRF 
AALI(J)  =  2XAAL(J) 


Form 
TEST=RFI-RF 


Punch 
EER.RF 


Punch 

FORMAT 

1601 


-^-Ofes)^ (Is   TEST  =  0> 


Set 
LNC=LNC+I 


Set   III  sO, 
AL  =  AAL(LNC) 
ALI  =  AAL  (LNC) 


9p_ 


Set 
LNC  = 


Is     ITESTI/RF 
YesPX  =  CRIT2 


Form 

TEST  =  ALI-AL 


Punch 
FORMAT    601 


-+ — presV-— (Is     TEST  =  6) 


Go    to    LOOP 
(See   next  page) 


Punch    SUM, 
SLOPE,  AL 


NO 


Is    ITESTI/ALX 

i    CRIT  I         A*iYes, 


Form      Augumented 
T    matrix    and 
Calculate     X(J) 


Calculate  TERM(LNC) 
BESEL(I,J),  COEF(I,J), 
T    matrix 


-(noV-(is    INDEX(LNC)  =  0> 


Punch    X(J) 


Calculate     ERR 


Punch    RF 
Punch    FLUX 
and    RADIUS 


f     SENSEX 
; WITCH    1/ 


I 

SWI 


Go    to    902 
(next    page) 


Calculate    f.0(  r.) 


Store     TLRM(LNC),  X(J) 
Set    AAL(LNC)=AL, 
AALI(LNC)=ALI 


Calculate    and 
Punch   Coefficients 


{ H0\+- <Is    LNC  =NC>— *-(jesV 


Set     INDEX(J)  =  0 
INDXI(J)  =  2 


91 


LOGIC     DIAGRAM    (CONT'D) 
APPENDIX       D 


ON 


If  MMM- 


Punch    RF,  ERR 
Set  EERR(3)  =  ERR 
RFI  =RRF(3)=RF 


LOOP 


If    SENSE 
SWITCH  4 


OFF 


Increment     RF 
Punch   new    RF 
Set    RRF(J)  =  RRF(J+|) 
EERR(J)  =  EERR(J+I) 


Punch   ERR,RF 
Set  RRF(I)  =  RF 
EERR(I)  =  ERR 


Punch    ERR,  RF 
Set  MMM  =  2 

RFI  =  RRF(2)=RF 
EERR(2)=ERR 
RF=RF  f  DLTRF 


Set    AALO(J)=AAL(J) 
Increment    AAA(J) 


Go  to    888 
(Preceding  page) 


'If    SENS 
.SWITCH  4 


ON 


Set   MMM=I,  RFI=RF 
RF  =  RF  +  DLTRF 

AALO(J)  =AAL(J) 
Increment    AAL(J) 


Go    To    80 
(Preceding  poge) 


Is     111  =  0 


©"* 


902 


Calculate 
SUM  =  ITI 


If    INDXKLNO-I 


Print     SUM 
Set  III  =  I 

ALI= AL 
SUMI =SUM 
AL=AL+DL 


Go    to    90 
(Preceding  page) 


NO 


^ 


Ki> 


Calculate 
SLOPE 


If    SENSE 
SWITCH   2 


ON 


Print    SUM 
SLOPE, AL 


Set  INDXI(LNC)=| 
SUMI=SUM,  ALI=AL 
AL=AL-SUM/SLOPE 
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C    LEAST  SQUARES  PROGRAM  USING  UP  TO  4  HARMONICS,  P3  APPROXIMATION 

C 

DIMENSION  T( 12,13) »X (12) , ROOT (8) ,TERM(8) ,ARG(20) 
DIMENSION  PHI(50),R(5C),W(50),ALSIG(4)»BTSIG(4)»ENSIG(4) 
DIMENSION  BESEL(50,4) ,GAM0(2) ,GAM1 (2),COEF(8»12),CFR(6) 
DIMENSION  OMSIGU)  »Z(36)  ,  TEEM  (32)  »AAL(4)  ,AALK4)  »INDEX(4) 
DIMENSION  SLOPE (4) ,AAL0(4) iRRF(3) ,EERR(3) ,INDX1(4) 

1  FORMAT(Elb.7,E16.7,E16.7,E16.7) 

2  FORMAT (E16. 7/) 

3  F0RMAT(E16.7,E16.7,E16.7/) 

4  FORMAT ( 15, 15,15, 15,15,15) 

5  F0RMAT(E16.7,E16.7,E16.7,E16.7/ ) 

6  FORMAT! E16. 7 , E 16. 7/ ) 

7  FORMAT(E14.7,I6,E14.7) 
C 

8  READ  4,NM0D»MMM 
PRINT  4,NM0D,MMM 
PUNCH  4,NM0D,MMM 

READ  l,GAMO(  1)  ,GAM0(2)  ,GAM1(  1  )  .GAMK2) 

PUNCH  1,GAM0( 1) ,GAM0(2) ,GAM1( 1 ) ,GAM1(2) 

READ  1,SIGF,SIGM,RM,CRIT1 

PUNCH  1,SIGF,SIGM,RM,CRIT1 

READ  1,CRIT2,FACTR,DIV1,DLTRF 

PUNCH  1  ,CRIT2,FACTR,DIV1,DLTRF 

READ   l.EERR(l) ,EERR(2) ,RRF ( 1 ) , RRF ( 2 ) 

PUNCH  l.EERR(l)  ,  EERR ( 2 )  ,RRF ( 1 ) , RRF ( 2 ) 

DO  15  J=l,NMOD 

READ  7, PHI (J) »K»W( J) 

15  W( J)=1./(W( J)**2> 
DO  16  J=l,NMOD 
READ  l.A 

16  R(J)=A*FACTR/DIV1 
READ  4,NC 

PUNCH  4»NC 
888  READ  1,RF,ADD1,DL 
DO  3115  J=1,NC 

READ  1,AAL( J) ,AAL0( J) ,SLOPE( J) 
READ  4,INDEX( J) 
AAL1(J)=2.*AAL( J) 
INDXK  J)=0 
3115  PUNCH  1 »AAL( J)  ,AALO( J)  ,SLOPE( J) 
PRINT  3»RF,ADD1,DL 
PUNCH  3»RF,ADDl,DL 
C 

J  =  l 

K  =  l 

12  Kl-J+1 

A  =  l.  +GAMU(K)*(.8  +  27.*GAMHK)/35.  ) 

ROOT(J)-  3  5.*A*(1.-SQRT(1.-108.*GAM0(K)*GAM1(K)/(35.*A**2)  )  )/18. 

RGOTf J+1)=35.*A* (l.  +  SQRTl 1 . -  1 08 . *GAM0 ( K ) *GAMl ( K ) / ( 35.*A**2)  ) )/18. 

ROOT(  J  +  2)  =35.*GAMKK)/(8.*GAM1  (K)+7.  ) 

ROOT( J+3)  =7. 

DO  9  1  = J  » K 1 
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V    CFRU  )  =  (ROOT(  1  l-GAMO  ( K  )  *  (  3  •*GAM1  (  K,  )  +  V.  )  )/(ROOT(  I  )*(ROCT(  I  )-7.  )  ) 
IF ( J-l Jo.lL. 11 

10  J  =  3 

K  =  2 

GO    TO    12 

11  RF1=2.*RF 

80  TEST=RF1~RF 
IF(TEST)81»82»83 

81  TEST=-TEST 

8  3  IF(TES1 /RF-CRIT2) 1316.900,900 

82  PRINT  1601 

1601  F0RMATU7H  DELTA  RADIUS  =  0/) 
GO  TO  1316 
C 

C      FROM  HERE  TO  1800  FINDS  EIGEN-VALUE  AND  E IGEN-VECTORS 
C 

900  DO  1800  LNC=1»NC 

111=0 

KKK  =  0 

AL=AAL(LNC) 

AL1=AAL1(LNC) 

90  TEST=AL1-AL 
IFITEST )91 ,92,93 

91  TEST=-TEST 

93  IF(TEST/AL-CRIT1)316.94.94 

92  PUNCH  601 

601  FORMAT* 17H  DELTA  LAMBDA  =  0/) 
GO  TO  316 

94  SL=AL**2 

110  DO  111  J=l,8 

IF( J-5 ) 112  .1  11  .  112 

112  TERM( J)=SORT(ROOT( JJ-SL) 

111  CONTINUE 
TERM(5)=SQRT(SL-ROOT(5) ) 
DO  113  J=l»4 

ARG( J)=TERM( J)*SIGF*RF 
ARG( J+4)=TERM( J+4)*SIGM*RF 
ARGIJ+8 )=ARG( J  +  4) 
ARG( J< 12) =TERM( J+4) *SIGM*RM 

113  ARG( J+16) -ARG( J+12) 
DO  114  J=l,20 
IF(J-5)115,118.119 

115  BESEL(J,1)-I2ER(ARG(J)  ) 

BESEL(J,2)  =  I0NE(ARG( J)  ) 

A-  -1. 
124  B-  -A 

GO  TO  123 

118  BE5EL(J»1)=YZER(ARG( J)  ) 
BLSEL( J,2)=Y0NE(ARG(J)  ) 
GO  TO  120 

119  IF(J-9) 121,116,122 

121  BESEL(J.1)=K2ER(ARG( J) ) 
BCGELl J.2 J=KONE(ARG( J)  ) 
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A  =  l. 

B  =  A 

CC  TO  123 

116  BE5EL(J.1)=JZER(ARG( J) ) 
BESEL(J»2)=JONE(ARG(J) ) 
IFU-9)  120.125.12C 

125  BESELU.2 )  =  I ARG ( J ) **2 ) * ( 1 . -ARG ( J ) **2/ 12 •  )/8. 

BE5ELU  »4>  =  (  ARG(  J)**3)*(  1 .  -ARG  (  J  )  **2/  16  .  )/48. 

GO  TO  114 
120  A=l. 

GO  TO  124 

122  IFU-13)  115,118,117 

117  I F ( J-17 ) 12  1  .  116  ,  115 

123  BE5ELI J ,3 ) =2 .*A*EESEL ( J, 2) /ARG( J)   +B*BESEL ( J  ,  1  ) 
BESEL(J»4)=4.*A*BESEL(J,3)/ARG(J)   +B*BESEL ( J  ,  2  ) 

114  CONTINUE 
A  =  -l. 
B  =  A 
C^-A 
L  =  l 
M  =  2 
K  =  L 

133  DO  130  J^L.M 
C0EF(1.J)=A 

C0EF(2,J) =B*6.*TERM( J)*AL*CFR(J) 

COEF(3»J)=A*3.*GAM0(K)*AL/ROOT( J) 

COEF (4,  J)=A*(3.*SL-R00T(  J)  )*CFR( J) 

COEF (5. J) =C*TERM( J ) *COEF ( 3 » J ) /AL 

C0EF(6.J) --B*?.*(SL-ROOT( J) )*CFR( J) 

C0EF(7,J) =8*9.*TERM(J)*(5.*SL-ROOT< J) )*CFR( J)/10. 

130  C0EF(8,J) -C*TERM< J)*C0EF(6. J) /2. 
IFIM-2) 132,132,131 

132  L  =  5 

K  =  2 

A  =  -A 
136  C=-C 

M  =  L 

GO  TO  133 

131  IF(M-5) 135,135,134 
135  L  =  6 

B  =  -B 

GO  TO  136 

134  K  =  l 

A  =  -l. 
J^3 
142  COEF (1,J) -0. 

COEF(2»  J)  =-A*5.*GAMl  (  K )  *  (  2  .*SL-ROOT  (  J  )  )/(  TERM(  J)*ROOT(  J)  ) 

COEF(3,J) =A 

COEF (4, J) =A*5.*GAMl(K)*AL/ROOT(J) 

COEF(5,J) =-A*AL/TERM(J) 

COEF(6,J) =-C0EF(4, J) 

COEF (7, J) =-A*GAMKK)*AL*( 15. *SL-1 1 .*ROOT ( J ) )/( 2#*ROOT ( J )*TERM( J ) ) 

COEF (8, J) =-TERM(J)*C0EF(4, J) /2. 
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CCEF(  l.J-f  1 
C0EF(2»J+1 
C0EF(3,J+1 
C0EF(4,J+1 
C0EF<5,J+1 
C0EF(6, J+l 
C0EF(7.Jt  -1 
C0EF(8,J+1 
IF( J-3) 14 1 

140  J=7 
A  =  -A 

GO  TO  142 

141  DO  146  K=5» 
DO  147  J=l, 

147  C0EF(J,k+4) 
IF(K-5) 146, 

148  DO  149  J=2, 

149  COEF(J,K)=- 
C0EF(7,K)=- 

146  CONTINUE 

DO  150  J=l, 
DO  151  K=l» 

151  T(K,J)=COEr 
T (4»J)=C0EF 
DO  152  K=2i 

152  T(K.J)=COEF 
Tt6,J)--C0EF 
T(7,J)=C0EF 

150  T(8,J)-^C0EF 
DO  155  J=9, 
DO  155  K=l , 

15  5  T ( J  ,  K ) - 0 . 0 
DC  156  J  =  5, 
T(9,J)  =COE 
T(9,J)  -  T( 
Tl  10,J)=2.« 
T (11,J)=-C0 
T( 11»J)=T( 1 
T( 12»J)=-C0 
T(12»J)=T(1 

156  T( 12»J)=T(  1 
IF( INDEX(LN 

902  5UM=CRAM(-1 
IF(  INDXKLN 

298  IF(II1 1300, 

300  111=1 

PRINT  602.  SUM 

602  FORMAT ( E 16.7 »6H 
AL1=AL 
SUM1=SUM 
AL=AL+DL 
GO  TO  90 


=  0.0 

=-A*2.*AL/TERM( J+l ) 

-  0 . 0 

=  A 

-0.0 

=-A*(SL+ROOT( J+l) )/(SL-ROOT( J+l) ) 

=  -A*( 3.*SL-ROOT( J+l )  )/(2.*TERM( J  +  l)  ) 

=A*(SL+ROOT( J+l ) )/ (2.* TERM (J+l ) ) 

140*141 


8 

8 

=COEF( J,K) 

146, 148 

8,3 

COEF( J.K ) 

COEF(7,K) 


12 
3,2 
(  K, 
(4, 
5,3 
(K, 
(  6, 
(  7, 


J)*BESEL( J»l ) 
J)*BESEL( J»l ) 

J)*BESELU»2  ) 
J)*BESEL(J»3 ) 
J)*BESEL( J»2 ) 

J)*BESEL( J, 4 ) 


1  2 
F(l 
9, J 
CO! 
EF( 
ltJ 
EF( 
2, J 

2  ,J 
C)  ) 
2.) 
O- 
300 


•J)*BESEL( J  +  8,1 )  +2.*COEF(5»J)*BE£EL( J+8.2J/3. 
)-(C0EF(4, j)*BESEL( J  +  0,1)  -COEF  (  6  ,  J  )  *BESEL  ( J  +  8,3)  )/8. 
F(3,J)*BESEL( J+8,1) /3.  +COEF ( 2 » J ) *6ESEL( J+8 »2 ) M. 
6»J)*BESEL/(J  +  8,3)M8.  +B.*COEF ( 7 , J ) *BESEL ( J+8 , 2 ) /2 1 . 
)  +(COEF(l,J)/4.  +7.*C0EF(4,J)/16. )*BESEL( J+8,1) 
liJ)*BESEL( J+8, 1 )/4.  +8.*COEF(8,J)*BESEL( J+8,4)/35. 
)  +COEF(A,J)*BFSEL( J+8,1 )/16. 
)  +3.*COEF<6, J)*BESEL( J+3,3)/16. 
902,902,352 

1  )298,?04,313 

,304 


SUM/  ) 


96 


3  04  SLOPE (LNC)  =  ( SUM-SUM  1 )/(AL-ALl ) 

IF (SENSE  SWITCH  2)330,313 
33C  IF(KKK)31*>314»31? 
3H  KKK  =  1 

PRINT  6  0-'+ 
604  FORMAT(8X4H  SUM.9X6H  SLOPE.11X7H  LAMbDA/) 

315  PRINT  3»SUM.SLOPE(LNC) »AL 
313  SUM 1= SUM 

INDX1(LNC)=1 
ALl^AL 

AL=AL-SUM/SLOPE(LNC) 
GO  TO  90 

316  PUNCH  6  0't 

PUNCH  3  .SUM. SLOPE (LNC) ,AL1 
352  DO  317  J= 1  *  1 2 

317  T( J, 13)=  -T( J.l ) 
DO  318  J=l»ll 

DO  318  1  =  1  .11 

318  T( J.I )=T( J.I+1 ) 
DO  319  J=l .11 
T(J»12)=0. 

319  T( 12»J)=0. 

X(l  )=CRAM( 12. ) 

IFISENSE  SWITCH  1)932.933 

932  PUNCH  610 

610  FORMAK5X7H  R(  1+1  )/) 

DO  935  J=l ,10 
935  PUNCH  1 »X( J) 

PUNCH  2.X  (  11  ) 

933  DO  1780  J=l,9 
K8=J  +8*(LNC-1) 
K9=J+  9*(LNC-1) 

IF( J-8) 1781 .1781.1 780 
1781  TEEM(K8)=TERM< J) 
178  0  Z(K9)~X(J) 

AAL1 (LNC) =AL1 
1800  AAL(LNC)~AL 

COMPUTES  LEAST  SQUARES  ERROR  AND  STARTS  ITERATION  ON  RF 

DO  2101  K=1»NC 

INDEX(K)-0 

INDX1 (K)=2 

J=l  +3*<K-1) 

ALSIGtK )=TEEM(J)*SIGF 

BTSIG<K)=TEEM( J+1)*SIGF 

ENSIG(K)=TEEM( J+4)*SIGM 
2  101  OMSIG(.".)-TFCM(  J+5)*SICM 

DO  235  K=l  ,NMOD 

IF(R(K)-RF)235. 235.236 
23  5  CONTINUE 
236  NROD=K-l 

DO  2001  J=l,12 
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20C1 


DO 

20G1 

T  (>. 

,K)  = 

DO 

240 

DO 

240 

L  =  l 

+  9* 

K  = 


1=1 ,NMOD 
(K-l) 
IF ( I-NROD) 241*241 #242 

241  BESELl I »K)=IZER(ALSIG(K)»R( I ) )  +  Z ( L) *IZER ( BTSIG ( K ) *R ( I ) ) 

GO  TO  240 

242  BESELIIiK) ~Z ( L  +  3 ) *YZER ( ENS  I G ( K ) *R ( I ) ) +Z ( L  +  4 ) *KZER ( OMS IG ( K ) *R ( I ) ) 
6ESCL( I »K)=BESEL( I »K)  +  Z  (  L  +  7  )  *  JZER  (  ENS  I G  (  K. )  *R  ( I ) ) 

BESEH I »K)=BESEL( I >K)  +Z ( L+8 ) * IZER (OMSIG ( K ) *R ( I ) ) 
24U  T(K*13)=T(K»13)  +  W ( I  ) *BESEL ( I »K  )  *PH I ( I ) 
DO  244  K=1.NC 
DO  244  J=K»NC 
DO  2  44  I=l,NMCu 

244  T(K»J)=T(K»J)  +W( I ) *BESEL ( I » J ) *BESEL( I ,K ) 
IF(NC-1 )8, 245, 246 

245  X( ]  )=T(1»13)/T(1»1) 
GO  TO  2248 

246  K2=NC-1 

DO  247  K=l *K2 

K1=K+1 

DO  247  J=K1.NC 

247  T(J»K)=T(KtJ) 
NC1=NC+1 

DO  248  J=NC1»12 

248  T(J,J)=1. 

X(l  )=CRAM( 12. ) 

2248  PUNCH  668 

668  FORMAT(2X12H  COEFICIENTS) 
DO  2249  J=1.NC 

2249  PUNCH  l.X(J) 

IFtSENSE  SWITCH  1)555.5555 
555  PUNCH  667, RF 
667  F0RMAT(E16.7,13H  =  ROD  RADIUS/) 

PUNCH  633 
633  FORMAT! 1X16H  CALCULATED  FLUX.5X7H  RADIUS/) 
5555  ERR=0. 

DO  2250  I=l,NMOD 
A  =  0. 

DO  2251  J=1,NC 
2251  A-A  +X( J)*PESEL( I »J) 

IF(SENSE  SWITCH  1)2520,2250 
2520  PUNCH  7, A, I ,R( I ) 

2250  ERR^EPR  +W( I ) * ( A-PHI ( I ) ) **2 
C 

C        SWITCH  4  OFF  FOR  ITERATION  ON  ROD  RADIUS 
C 

IFtSENSE  SWITCH  4)1200,1199 

1199  IFCMMM-1) 1 20U, 1201 t 1202 

1200  RRFf 1)=RF 
EERRt 1 J=ERR 
PUNCH  666, ERR 
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666  FCRV,AT(El6.7t9H  =  ERROR/) 

PUNCH  667, RF 

IF(5ENSE  SWITCH  4)888,1198 
119  8  MMM=1 

RF1=RF 

RF=RF  +DLIRF 

DO  1900  J=1»NC 

AAL0(J)=AAL( J) 
19C0  AAL(J)=AAL(J)*(1.  +ADD1 ) 

GO  TO  30 

1201  MMM=2 
RF1=RF 
RRF(2)=RF 
EERR(2)-ERR 
PUNCH  666, ERR 
PUNCH  667, RF 
RF=RF  +DLTRF 

1119  DO  1901  J=1,NC 

SLP=(AAL( J)-AAL0( J)  ) / ( RRF( 2 ) -RRF ( 1)  ) 

AAL0( J)=AAL( J) 
1901  AAL(J)=AAL(J)  +  (  RF-RRF  (  2  )  )  *SI_P 

GO  TO  80 

1202  RRF(3)=RF 
RF1=RF 
EERR(3)=ERR 
PUNCH  667, RF 
PUNCH  666, ERR 
DO  1907  J-1,12 
DO  1907  K=l»13 

1907  T(J,K)=0, 

DO  1908  J=4.12 

1908  T(J,J)=1. 

DC  1910  J=l,3 
T( J.l  )=1. 
T(J,2)=RRF( J) 
T( J,3)=RRF( J)**2 

1910  T( J,13)=EERR< J) 
X( 1 )=CRAM( 12.  ) 

RF=  -X(2)/(2.*X(3)  ) 
PUNCH  667, RF 
DO  1911  J  =  l,2 
EERR(J)=EERR( J+l) 

1911  RRF(J)=F,RF(  J  +  l  ) 
GC  TO  1119 

1316  PUNCH  666,ERR 
PUNCH  667, RF 
GO  TO  8 
END 
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A  study  Is  made  of  the  possibility  of  using  a  fictional 
"equivalent  rod  radius"  to  accurately  predict  the  neutron  flux 
depression  in  strong  absorbers  by  the  use  of  the  P^  approx- 
imation to  the  Boltzmann  equation  for  monoenergetic  neutron 
transport. 

The  theory,  appropriate  boundary  conditions,  and  the  re- 
quired computer  programs  were  developed  for  determining  the 
equivalent  rod  radius  from  experimental  measurements  in  assem-  - 
blies  having  exponential  flux  dependence  in  the  z-  direction. 
Only  cylindrical  geometry  is  considered  in  the  theoretical 
development. 

Application  of  the  commonly  applied  boundary  conditions 
yielded  16  equations  with  which  to  determine  12  unknowns.  A 
study  of  possible  boundary  conditions  was  made,  and  one  com- 
bination of  12  equations  was  found  to  give  a  good  approximation 
to  the  desired  physical  conditions.  It  is  felt  that  the  same 
seleotion  rule  can  be  used  to  determine  "correct"  boundary 
conditions  for  higher  order  PL  approximations,  when  L  is  odd. 

Although  much  of  the  experimental  work  is  still  to  be 
accomplished,  the  preliminary  experimental  work  performed  here 
strongly  indicates  that  the  determination  of  good  "equivalent 
radii"  is  entirely  feasible. 


